approximation by Lax, see (4.2), when the assumption (4.2) is relaxed:
$\mbox{$\langle$}\overline{\mbox{\boldmath$\zeta$}}_{k}\mbox{$\rangle$}(\mbox{\boldmath$\bf
z$})\not\equiv{\rm const}$ at $\mbox{\boldmath$\bf z$}\in v_{k}^{l},$
$(k=i,j$; a case of the local elasticity was considered in [44], [45]) that is
equivalent to the equality

$\displaystyle\mbox{$\langle$}\mbox{\boldmath$\cal
L$}^{\theta\zeta}_{j}(\mbox{\boldmath$\bf
z$},\overline{\mbox{\boldmath$\zeta$}}_{j})|;v_{1},{\bf
x}_{1}\mbox{$\rangle$}$ $\displaystyle=$
$\displaystyle\mbox{$\langle$}\mbox{\boldmath$\cal
L$}^{\theta\zeta}_{j}(\mbox{\boldmath$\bf
z$},\overline{\mbox{\boldmath$\zeta$}}_{j})\mbox{$\rangle$}.$ (4.12)

It greatly simplifies the problems (4.5) and (4.1) where in such a case
$\mbox{\boldmath$\cal J$}^{\theta\zeta}_{i,j}=\mbox{\boldmath$\cal
J$}_{i,j}^{\theta\zeta\infty}$ ($\mbox{\boldmath$\bf x$}\in v_{i}$) reducing
these equations to

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\langle\overline{\mbox{\boldmath$\vartheta$}}_{i}^{[n+1]}\rangle(\mbox{\boldmath$\bf
z$})\\!$ $\displaystyle=$
$\displaystyle\\!\langle\mbox{\boldmath$\vartheta$}\rangle({\bf
z})\\!+\\!\\!\int\\!\\!\\!\mbox{\boldmath$\cal
J$}^{\theta\zeta\infty}_{i,j}(\mbox{$\langle$}\overline{\mbox{\boldmath$\zeta$}}_{j}^{[n]}\mbox{$\rangle$})(\mbox{\boldmath$\bf
z$})[\varphi(v_{j},{\bf x}_{j}|;v_{i},{\bf x}_{i})-n^{(j)}(\mbox{\boldmath$\bf
x$}_{j})]d{\bf x}_{j},$ (4.13)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\langle\overline{\mbox{\boldmath$\vartheta$}}_{i}^{[n+1]}\rangle(\mbox{\boldmath$\bf
z$})\\!$ $\displaystyle=$
$\displaystyle\\!\mbox{\boldmath$\vartheta$}^{b(0)}({\bf
z})+\\!\\!\int\\!\\!\mbox{\boldmath$\cal
J$}^{\theta\zeta\infty}_{i,j}(\mbox{$\langle$}\overline{\mbox{\boldmath$\zeta$}}_{j}^{[n]}\mbox{$\rangle$})(\mbox{\boldmath$\bf
z$})\varphi(v_{j},{\bf x}_{j}|;v_{i},{\bf x}_{i})\,d{\bf x}_{j},$ (4.14)

respectively. For local elasticity, this method uses as a background the new
GIE proposed in [44], [45] and permits to abandonment of the basic concepts of
micromechanics: the hypothesis of “ellipsoidal symmetry”, and the effective
field hypothesis (EFH H1a considered in the next Section 5.2). Some new
effects were discovered that were impossible in the the framework of the
classical (the first) background of micromechanics.

Equations (4.9) and (4.10) are obtained for general nonlinear cases of either
the state-based (2.8) or bond-based (e.g. (2.15)) PM. Equation (4.13) and
(4.14) are reduced to the corresponding equations obtained before for either
linear bond-based (2.10) PM [57], [58] or linear state-based PM (2.13) and
(2.14) [64].

###  4.2 Body forces with compact support

The scheme (4.9) forms the Neumann series for the solution
($\mbox{\boldmath$\bf x$}\in v_{i}^{l}$)

$\langle\overline{\mbox{\boldmath$\bf
u$}}\mbox{$\rangle$}_{i}(\mbox{\boldmath$\bf
x$}):=\langle\lim_{n\to\infty}\overline{\mbox{\boldmath$\bf
u$}}^{[n+1]}\mbox{$\rangle$}_{i}(\mbox{\boldmath$\bf
x$})=\widehat{\mbox{\boldmath$\bf D$}}_{i}^{ub}(\mbox{\boldmath$\bf
b$},\mbox{\boldmath$\bf x$}),$ (4.15)

which yields the representation for the statistical average of displacements
inside the fixed extended inclusion $(\mbox{\boldmath$\bf x$}\in v_{i}^{l})$
(conditional average)

$\mbox{$\langle$}{\mbox{\boldmath$\bf
u$}}\mbox{$\rangle$}_{i}(\mbox{\boldmath$\bf x$})=\widehat{\mbox{\boldmath$\bf
D$}}_{i}^{ub}(\mbox{\boldmath$\bf b$},\mbox{\boldmath$\bf
x$})+\mbox{\boldmath$\cal L$}_{i}^{uu}(\mbox{\boldmath$\bf
x$}-\mbox{\boldmath$\bf x$}_{i},\widehat{\mbox{\boldmath$\bf
D$}}_{i}^{ub}(\mbox{\boldmath$\bf b$},\mbox{\boldmath$\bf x$})).$ (4.16)

It should be mentioned that the tensor $\widehat{\mbox{\boldmath$\bf
D$}}_{q}^{ub}(\mbox{\boldmath$\bf b$},\mbox{\boldmath$\bf y$})$ (4.15) is an
inhomogeneous function of coordinates of the moving inclusion
$\mbox{\boldmath$\bf y$}\in v_{q}^{l}$. Moreover, the tensors
$\widehat{\mbox{\boldmath$\bf D$}}_{q}^{ub}(\mbox{\boldmath$\bf
b$},\mbox{\boldmath$\bf y$})$ (4.15) depend on all interecting inclusions (at
least at $|\mbox{\boldmath$\bf x$}_{q}|\leq a^{\delta}+l_{\delta}$).

Equation (4.16) allows us to estimate a statistical average of the
displacement field in the macropoint $\bf X$

$\langle{\mbox{\boldmath$\bf u$}}\mbox{$\rangle$}(\mbox{\boldmath$\bf
X$}):=c^{l(0)}\mbox{$\langle$}\mbox{\boldmath$\bf
u$}\mbox{$\rangle$}^{l(0)}(\mbox{\boldmath$\bf
X$})+c^{l(1)}\mbox{$\langle$}\mbox{\boldmath$\bf
u$}\mbox{$\rangle$}^{l(1)}(\mbox{\boldmath$\bf X$})$ (4.17)

expressed through the statistical averages of displacements in the point $\bf
X$ in the matrix $\mbox{$\langle$}\mbox{\boldmath$\bf
u$}\mbox{$\rangle$}^{l(0)}(\mbox{\boldmath$\bf X$})$ and inclusions
$\mbox{$\langle$}\mbox{\boldmath$\bf
u$}\mbox{$\rangle$}^{l(1)}(\mbox{\boldmath$\bf X$})$. Indeed, at first we
built some auxiliary set $v^{1}_{i}(\mbox{\boldmath$\bf X$})$ with the
indicator function $V^{1}_{i}(\mbox{\boldmath$\bf X$})$ and the boundary
$\partial v^{1}_{i}(\mbox{\boldmath$\bf X$})$ formed by the centers of
translated ellipsoids $v_{q}({\bf 0})$ around the fixed point $\bf X$. We
construct $v^{1}_{i}(\mbox{\boldmath$\bf X$})$ as a limit $v^{0}_{ki}\to
v_{q}^{1}(\mbox{\boldmath$\bf X$})$ if a fixed ellipsoid $v_{k}$ is shrinking
to the point $\bf X$. Then $\mbox{$\langle$}\mbox{\boldmath$\bf
u$}\mbox{$\rangle$}^{l(1)}(\mbox{\boldmath$\bf X$})$ can be estimated as
($\mbox{\boldmath$\bf y$}\in v^{l}_{q})$

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!c^{l(1)}\mbox{$\langle$}\mbox{\boldmath$\bf
u$}\mbox{$\rangle$}^{l(1)}(\mbox{\boldmath$\bf X$})\\!$ $\displaystyle=$
$\displaystyle\\!c^{l(1)}\mbox{\boldmath$\bf u$}^{b(0)}(\mbox{\boldmath$\bf
X$})+\\!\\!\int_{v^{1}_{i}({\bf
X})}\\!\\!\\!\\!\\!\\!\\!n^{(1)}(\mbox{\boldmath$\bf
x$}_{q})\mbox{\boldmath$\cal L$}^{uu}_{q}(\mbox{\boldmath$\bf
X$}\\!-\mbox{\boldmath$\bf x$}_{q},\mbox{\boldmath$\cal
D$}_{q}^{u}(\mbox{\boldmath$\bf b$},\mbox{\boldmath$\bf y$}))d{\bf x}_{q},$
(4.18)

and, therefore, the statistical average of displacements (4.17) is expressed
through the body force density ($\mbox{\boldmath$\bf y$}\in v^{l}_{q})$

$\displaystyle\mbox{$\langle$}\mbox{\boldmath$\bf
u$}\mbox{$\rangle$}(\mbox{\boldmath$\bf X$})\\!$ $\displaystyle=$
$\displaystyle\\!\mbox{\boldmath$\bf u$}^{b(0)}(\mbox{\boldmath$\bf
X$})+c^{(0)}\int\mbox{\boldmath$\cal L$}^{uu}_{q}(\mbox{\boldmath$\bf
X$}-\mbox{\boldmath$\bf x$}_{q},\mbox{\boldmath$\cal
D$}_{q}^{ub}(\mbox{\boldmath$\bf b$},\mbox{\boldmath$\bf
y$}))\varphi(v_{q},{\bf x}_{q}|;\mbox{\boldmath$\bf X$})d{\bf x}_{q}$ (4.19)
$\displaystyle+$ $\displaystyle\\!\int_{v^{1}_{i}({\bf
X})}n^{(1)}(\mbox{\boldmath$\bf x$}_{q})\mbox{\boldmath$\cal
L$}^{uu}_{q}(\mbox{\boldmath$\bf X$}-\mbox{\boldmath$\bf
x$}_{q},\mbox{\boldmath$\cal D$}_{q}^{ub}(\mbox{\boldmath$\bf
b$},\mbox{\boldmath$\bf y$}))d{\bf x}_{q},$

where the first and the second integrals in Eq. (4.18) correspond to the first
and the second terms in the right-hand side of Eq. (4.17), and
$\varphi(v_{q},{\bf x}_{q}|;\mbox{\boldmath$\bf X$})=0$ if
$\mbox{\boldmath$\bf x$}_{q}\in v^{1}_{i}(\mbox{\boldmath$\bf X$})$.

The macroscopic stresses in the effective constitutive law are

$\mbox{$\langle$}\mbox{\boldmath$\sigma$}\mbox{$\rangle$}(\mbox{\boldmath$\bf
x$})=\mbox{$\langle$}\mbox{\boldmath$\cal L$}^{\sigma(0)}(\mbox{\boldmath$\bf
u$})\mbox{$\rangle$}(\mbox{\boldmath$\bf
x$})+\mbox{$\langle$}\mbox{\boldmath$\tau$}\mbox{$\rangle$}(\mbox{\boldmath$\bf
x$})$ (4.20)

where the first item on the right-hand side is simplified for the linear
matrix (see Eq. (3.43)) and depends on the statistical averages
$\mbox{$\langle$}\mbox{\boldmath$\bf u$}\mbox{$\rangle$}(\mbox{\boldmath$\bf
x$})$ ($\mbox{\boldmath$\bf x$}\in\mathbb{R}^{d}$) (4.19). For an estimation
of an average
$\mbox{$\langle$}\mbox{\boldmath$\tau$}\mbox{$\rangle$}(\mbox{\boldmath$\bf
x$})$, we consider, at first, the fixed inclusion $v_{q}$ with the center
$\mbox{\boldmath$\bf x$}_{q}$. This inclusion produces value
$\mbox{\boldmath$\cal L$}^{\sigma}_{1}(\mbox{\boldmath$\bf
x$}-\mbox{\boldmath$\bf x$}_{q},\mbox{\boldmath$\bf u$})=\mbox{\boldmath$\cal
L$}^{\sigma}(\widetilde{\mbox{\boldmath$\tau$}})$ estimated by Eq. (2.17) with
replacement $\widetilde{}\mbox{\boldmath$\bf
f$}^{(0)}\to\widetilde{}\mbox{\boldmath$\bf f$}_{1}$ as in Eq. (3.26). Then a
statistical average of the local polarization tensor
$\mbox{$\langle$}\mbox{\boldmath$\tau$}\mbox{$\rangle$}(\mbox{\boldmath$\bf
X$})$ is also obtained by averaging over the domain $v_{i}^{1}$

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\mbox{$\langle$}\mbox{\boldmath$\tau$}\mbox{$\rangle$}(\mbox{\boldmath$\bf
X$})$ $\displaystyle=$ $\displaystyle\int_{v^{1}_{i}({\bf
X})}n^{(1)}\mbox{\boldmath$\cal L$}^{\sigma}(\mbox{\boldmath$\bf
C$}_{1}\mbox{\boldmath$\eta$}^{D})(\mbox{\boldmath$\bf y$}-\mbox{\boldmath$\bf
X$})d\mbox{\boldmath$\bf y$},$ (4.21)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\mbox{\boldmath$\eta$}^{D}(\mbox{\boldmath$\bf
x$},\mbox{\boldmath$\bf y$})$ $\displaystyle=$
$\displaystyle\Big{[}\mbox{\boldmath$\cal D$}_{i}^{ub}(\mbox{\boldmath$\bf
b$},\mbox{\boldmath$\bf x$})+\mbox{\boldmath$\cal
L$}_{i}^{uu}(\mbox{\boldmath$\bf x$}-\mbox{\boldmath$\bf
x$}_{i},\mbox{\boldmath$\cal D$}_{i}^{ub}(\mbox{\boldmath$\bf
b$},\mbox{\boldmath$\bf x$})$ (4.22) $\displaystyle-$
$\displaystyle\mbox{\boldmath$\cal D$}_{i}^{ub}(\mbox{\boldmath$\bf
b$},\mbox{\boldmath$\bf y$})-\mbox{\boldmath$\cal
L$}_{i}^{uu}(\mbox{\boldmath$\bf y$}-\mbox{\boldmath$\bf
x$}_{i},\mbox{\boldmath$\cal D$}_{i}^{ub}(\mbox{\boldmath$\bf
b$},\mbox{\boldmath$\bf y$}))\Big{]},$

where $v_{q}^{1}(\mbox{\boldmath$\bf x$})$ is defined by (2.34) and one used
the representation for the statistical average of the displacement field in
the inclusion $\mbox{\boldmath$\bf x$}\in v_{q}$:
$\mbox{$\langle$}\mbox{\boldmath$\bf
u$}\mbox{$\rangle$}_{q}(\mbox{\boldmath$\bf x$})$ (4.18).

New nonlocal effective constitutive low (4.22) looks very similar to the
nonlocal effective constitutive law

$\mbox{$\langle$}\mbox{\boldmath$\sigma$}\mbox{$\rangle$}(\mbox{\boldmath$\bf
x$})={{}^{L}}\\!\mbox{\boldmath$\bf
L$}^{(0)}\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}(\mbox{\boldmath$\bf
x$})+\mbox{$\langle$}^{L}\\!\mbox{\boldmath$\tau$}\mbox{$\rangle$}(\mbox{\boldmath$\bf
y$}),\ \ \ ^{L}\\!\mbox{\boldmath$\tau$}(\mbox{\boldmath$\bf
y$}):={{}^{L}}\\!\mbox{\boldmath$\bf L$}_{1}(\mbox{\boldmath$\bf
y$})\mbox{\boldmath$\varepsilon$}(\mbox{\boldmath$\bf y$}),$ (4.23)

corresponding to locally elastic constitutive law (2.2) at the microlevel.
But, the first term ${{}^{L}}\\!\mbox{\boldmath$\bf
L$}^{(0)}\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}(\mbox{\boldmath$\bf
x$})$ in the right-hand side of Eq. (4.23) corresponds to the arbitrary
inhomogeneous
$\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}(\mbox{\boldmath$\bf
x$})\not\equiv$const. whereas $\mbox{$\langle$}\mbox{\boldmath$\cal
L$}^{\sigma(0)}(\mbox{\boldmath$\bf u$})\mbox{$\rangle$}(\mbox{\boldmath$\bf
x$})$ (4.23) is reduced to $\\!\mbox{\boldmath$\bf
L$}^{(0)}\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}(\mbox{\boldmath$\bf
x$})$ only for homogeneous
$\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}(\mbox{\boldmath$\bf
x$})\equiv$cons.

Numerical analysis was performed for 1D case [62], [63] for selfequilibrated
and vanished outside some loading region body force: $\mbox{\boldmath$\bf
b$}(\mbox{\boldmath$\bf x$})=0$ for $|\mbox{\boldmath$\bf x$}|>a^{b}$. The
engineering approach using the scale separation hypothesis (2.27) can be
described in the following scheme: estimation of $\mbox{\boldmath$\bf L$}^{*}$
(4.18) at the homogeneous VBC (2.31), and evaluation of
$\mbox{$\langle$}\mbox{\boldmath$\bf u$}^{\rm
EA}\mbox{$\rangle$}(\mbox{\boldmath$\bf x$})$ by Eq. (2.1)-(2.2) with
replacement $\mbox{\boldmath$\bf L$}\to\mbox{\boldmath$\bf L$}^{*}$. It leads
to monotonical $\mbox{$\langle$}\mbox{\boldmath$\bf u$}^{\rm
EA}\mbox{$\rangle$}(\mbox{\boldmath$\bf x$})$ whereas
$\mbox{$\langle$}\mbox{\boldmath$\bf u$}\mbox{$\rangle$}(\mbox{\boldmath$\bf
x$})$ (6.6) (simultaneous consideration of peridynamic constitutive equation
(2.4) and inhomogeneity of $\mbox{\boldmath$\bf b$}(\mbox{\boldmath$\bf x$})$)
is strongly nonmonotonical dirrering from $\mbox{$\langle$}\mbox{\boldmath$\bf
u$}^{\rm EA}\mbox{$\rangle$}(\mbox{\boldmath$\bf x$})$ by a factor 9.

For the linear bond-based model of constituents, Eqs. (4.22)-(4.22) are
reduced to the corresponding equations in [62], [63]. The representations
(4.23) and (4.22) can be generalized to the nonlinear inclusion properties
analogous to Eqs. (4.1). A nonlocal effective operator similar to (4.22) was
also obtained from Eq. (4.16) for the linear bond-based model in [62];
generalization to the linear state-based model is obvious.

Comment 4.1. It is interesting that to the best authors’ knowledge, there are
no academic publications considering the mentioned effective deformations of
locally elastic composites produced by prescribed inhomogeneous body forces
although the practical importance of this problem is well known, e.g., for
laser heating of CMs widely used in both the additive manufacturing [298] and
military applications [147]. So, laser pulse heating of surfaces results in
the rapid rise of temperature in the region irradiated by the laser beam. This
produces a high-temperature gradient with high thermal stress levels in this
region. A classical approach is based on the estimation of thermal stresses in
a homogeneous half-space subjected to laser heating (see, e.g., [300]). In so
doing, the implementation of the additive technologies of laser heating can
produce the local fields of both the temperature and thermal stresses with a
scale comparable to the material scale of particles in CMs [298] that
indicates the need to analyze the nonlocal effects (as in Eqs. (4.16) and
(4.19)).

###  4.3 Estimation of a set of surrogate operators

The loading by inhomogeneous body force $\mbox{\boldmath$\bf
b$}(\mbox{\boldmath$\bf x$})$ can also be exploited as a very prospective tool
for data-driven learning of the surrogate nonlocal constitutive laws of CMs
rather than just a particular parameter of loading. This perspective direction
initiated by Silling [258], [302] for periodic systems can be generalized to
the random structure CM in the following way. So, we define two datasets

$\displaystyle{\cal D}$ $\displaystyle=$
$\displaystyle\\{\mbox{$\langle$}\mbox{\boldmath$\sigma$}_{k}\mbox{$\rangle$}(\mbox{\boldmath$\bf
b$}_{k},\mbox{\boldmath$\bf x$}),\mbox{$\langle$}{\mbox{\boldmath$\bf
u$}}_{k}\mbox{$\rangle$}(\mbox{\boldmath$\bf b$}_{k},\mbox{\boldmath$\bf
x$}),\mbox{$\langle$}{\mbox{\boldmath$\bf
u$}}_{ik}\mbox{$\rangle$}(\mbox{\boldmath$\bf b$}_{k},\mbox{\boldmath$\bf
x$}),\mbox{\boldmath$\bf b$}_{k}(\mbox{\boldmath$\bf x$})\\}_{k=1}^{N},$
(4.24) $\displaystyle\widetilde{\cal D}$ $\displaystyle=$
$\displaystyle\\{\widetilde{\mbox{\boldmath$\sigma$}}_{k}(\mbox{\boldmath$\bf
b$}_{k},\mbox{\boldmath$\bf x$}),\widetilde{\mbox{\boldmath$\bf
u$}}_{k}(\mbox{\boldmath$\bf b$}_{k},\mbox{\boldmath$\bf
x$}),\widetilde{\mbox{\boldmath$\bf u$}}_{ik}(\mbox{\boldmath$\bf
b$}_{k},\mbox{\boldmath$\bf x$}),\mbox{\boldmath$\bf
b$}_{k}(\mbox{\boldmath$\bf x$})\\}_{k=1}^{N},$ (4.25)

where each effective parameters
$\mbox{$\langle$}{\mbox{\boldmath$\sigma$}}_{k}\mbox{$\rangle$}(\mbox{\boldmath$\bf
x$})$, $\mbox{$\langle$}{\mbox{\boldmath$\bf
u$}}_{k}\mbox{$\rangle$}(\mbox{\boldmath$\bf x$})$ and
$\mbox{$\langle$}\mbox{\boldmath$\bf
u$}_{ik}\mbox{$\rangle$}(\mbox{\boldmath$\bf
x$}):=\mbox{$\langle$}\mbox{\boldmath$\bf
u$}_{k}\mbox{$\rangle$}_{i}(\mbox{\boldmath$\bf x$})$ (4.24) are estimated for
each $\mbox{\boldmath$\bf b$}_{k}$ in the offline stage in $k$-th realization;
the surrogate dataset $\widetilde{\cal D}$ is close to the dataset ${\cal D}$
of the effective field parameters presenting, in fact, a combination of both
the generation and compression stages; dataset ${\cal D}$ is called field PM
dataset. Dataset $\cal D$ defined by statistical averages in the
micromechanical problem is approximated by a surrogate model

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\mbox{\boldmath$\cal
L$}_{\rm\gamma}[\mbox{$\langle$}{\mbox{\boldmath$\bf
u$}}_{k}\mbox{$\rangle$}](\mbox{\boldmath$\bf x$})\\!\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\mbox{\boldmath$\Gamma$}(\mbox{\boldmath$\bf x$}),$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\mbox{\boldmath$\cal
L$}_{\rm\gamma}[\mbox{$\langle$}{\mbox{\boldmath$\bf
u$}}_{k}\mbox{$\rangle$}](\mbox{\boldmath$\bf x$})\\!\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\int\mbox{\boldmath$\bf
K$}_{\gamma}(|\mbox{\boldmath$\bf x$}-\mbox{\boldmath$\bf
y$}|)(\mbox{$\langle$}{\mbox{\boldmath$\bf
u$}}_{k}\mbox{$\rangle$}(\mbox{\boldmath$\bf
y$})-\mbox{$\langle$}{\mbox{\boldmath$\bf
u$}}_{k}\mbox{$\rangle$}(\mbox{\boldmath$\bf x$}))~{}d\mbox{\boldmath$\bf
y$},$ (4.26)

where $\gamma:=b,{\rm u_{i}},{\rm\sigma}$ and
${\mbox{\boldmath$\Gamma$}}(\mbox{\boldmath$\bf x$}):=-\mbox{\boldmath$\bf
b$}(\mbox{\boldmath$\bf x$}),\mbox{$\langle$}\mbox{\boldmath$\bf
u$}\mbox{$\rangle$}_{i}(\mbox{\boldmath$\bf
x$}),\mbox{$\langle$}\mbox{\boldmath$\sigma$}\mbox{$\rangle$}(\mbox{\boldmath$\bf
x$})$, respectively, correspond to three surrogate operators
$\mbox{\boldmath$\cal L$}_{\gamma}$. Obtaining an optimal surrogate model for
the kernel functions $\mbox{\boldmath$\bf K$}^{*}_{\gamma}$ corresponding to
the surrogate dataset $\widetilde{\cal D}$ is achieved by solving three
optimization problems

$\displaystyle\\!\\!\\!\mbox{\boldmath$\bf K$}_{\gamma}^{*}={\rm
arg}\\!\min_{\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!{{\bf
K}_{\gamma}}}\\!\sum_{k=1}^{N}\\!||\mbox{\boldmath$\cal
L$}_{\rm{\gamma}}[\mbox{$\langle$}{\mbox{\boldmath$\bf
u$}_{k}}(\mbox{\boldmath$\bf b$}_{k})\mbox{$\rangle$}](\mbox{\boldmath$\bf
x$})-\mbox{\boldmath$\Gamma$}_{k}(\mbox{\boldmath$\bf x$})||^{2}_{l_{2}}+{\cal
R}(\mbox{\boldmath$\bf K$}_{\gamma}),$ (4.27)

by the use, e.g., Adam optimizer [166] after each step of gradient descent,
see for details [302] (see also [118], [304], and [303]) where the kernel $\bf
K$ is presented as a linear combination of Bernstein-based polynomials. Here,
the $l_{2}$ norm is taken over $\mbox{\boldmath$\bf x$}\in w$, and ${\cal
R}(\mbox{\boldmath$\bf K$}_{\gamma})$ is a regularization term on the
coefficients that improves the conditioning of the optimization problem (e.g.,
Tikhonov regularization).

It should be mentioned that Eqs. (4.16) and (4.19) obtained for general
nonlinear cases of either the state-based (2.8) or bond-based (e.g. (2.15)) PM
are reduced to their linear counterparts of bond-based PM [62], [63]. However,
construction of the surrogate model (4.26) was considered in [62] only for the
linear bond-based model.

Of course, the scheme (4.24)-(4.26) can be generalized to the periodic media
in the spirit of Eq. (3.53). We consider a parallel transition of the grid
$\mbox{\boldmath$\bf x$}_{\alpha}$ over the vector $\chi$ with uniform random
distribution on $\Omega_{00}$. Then for each
$\mbox{\boldmath$\chi$}\in\Omega_{00}$ and fixed $\mbox{\boldmath$\bf
b$}^{k}$, we can estimate the statistical averages

$\displaystyle\\!\\!\\!\\!\mbox{$\langle$}\mbox{\boldmath$\sigma$}\mbox{$\rangle$}(\mbox{\boldmath$\bf
b$}_{k},\mbox{\boldmath$\bf x$})\\!$ $\displaystyle=$
$\displaystyle\\!\frac{1}{|\Omega_{00}|}\int_{\Omega_{00}}\\!\\!\\!\\!\mbox{\boldmath$\sigma$}(\mbox{\boldmath$\bf
b$}_{k},\mbox{\boldmath$\bf
x$},\mbox{\boldmath$\chi$})~{}d\mbox{\boldmath$\chi$},\ \
\mbox{$\langle$}\mbox{\boldmath$\bf u$}\mbox{$\rangle$}(\mbox{\boldmath$\bf
b$}_{k},\mbox{\boldmath$\bf
x$})=\frac{1}{|\Omega_{00}|}\int_{\Omega_{00}}\\!\\!\\!\\!\mbox{\boldmath$\bf
u$}(\mbox{\boldmath$\bf b$}_{k},\mbox{\boldmath$\bf
x$},\mbox{\boldmath$\chi$})~{}d\mbox{\boldmath$\chi$},$
$\displaystyle\\!\\!\\!\\!\mbox{$\langle$}\mbox{\boldmath$\bf
u$}\mbox{$\rangle$}_{i}(\mbox{\boldmath$\bf b$}_{k},\mbox{\boldmath$\bf
x$})\\!$ $\displaystyle=$
$\displaystyle\\!\frac{1}{|\Omega_{00}|}\int_{\Omega_{00}}\mbox{\boldmath$\bf
u$}_{i}(\mbox{\boldmath$\bf b$}_{k},\mbox{\boldmath$\bf
x$},\mbox{\boldmath$\chi$})~{}d\mbox{\boldmath$\chi$},$ (4.28)

where $\mbox{\boldmath$\sigma$}(\mbox{\boldmath$\bf
b$}_{k},\mbox{\boldmath$\bf x$},\mbox{\boldmath$\chi$}),\ \mbox{\boldmath$\bf
u$}(\mbox{\boldmath$\bf b$}_{k},\mbox{\boldmath$\bf
x$},\mbox{\boldmath$\chi$})$ and $\mbox{\boldmath$\bf
u$}_{i}(\mbox{\boldmath$\bf b$}_{k},\mbox{\boldmath$\bf
x$},\mbox{\boldmath$\chi$})$ are the DNS solutions obtained for each
$\mbox{\boldmath$\bf b$}_{k}(\mbox{\boldmath$\bf x$}),\mbox{\boldmath$\chi$}$.
The fields (4.28) are used for generation of dataset $\cal D$ (4.24). The rest
Eqs. (4.26) and (4.27) are not altered. It should be mentioned that the
authors [296] (see also [297]) and [304] used the averaged displacement field
$\mbox{$\langle$}\widetilde{\bf
u}\mbox{$\rangle$}^{\Omega}_{i}(\mbox{\boldmath$\bf
x$})=|\Omega_{i}|^{-1}\int_{\Omega_{i}}\widetilde{}\mbox{\boldmath$\bf
u$}(\mbox{\boldmath$\chi$})~{}d\mbox{\boldmath$\chi$}$ for the fixed grid
$\mbox{\boldmath$\bf x$}_{i}$ that leads to obtaining a discrete kernel
$\mbox{\boldmath$\bf K$}(|\mbox{\boldmath$\bf x$}|)$ as opposite to the
continuous ones obtained by averaging over the “moving averaging” cell (4.28)
(see for details p. 895 in [60]).

So, the offline dataset ${\cal D}$ (4.24) (and (4.24) with (4.28)) are formed
offline from the statistical averages $\mbox{$\langle$}{\mbox{\boldmath$\bf
u$}}_{k}\mbox{$\rangle$}(\mbox{\boldmath$\bf x$})$,
$\mbox{$\langle$}\mbox{\boldmath$\bf
u$}_{ik}\mbox{$\rangle$}(\mbox{\boldmath$\bf
x$}):=\mbox{$\langle$}\mbox{\boldmath$\bf
u$}_{k}\mbox{$\rangle$}_{i}(\mbox{\boldmath$\bf x$})$ and
$\mbox{$\langle$}{\mbox{\boldmath$\sigma$}}_{k}\mbox{$\rangle$}(\mbox{\boldmath$\bf
x$})$ (4.28) presenting, in fact, a combination of both the generation and
compression stages. Loading by dynamic body force $\mbox{\boldmath$\bf
b$}(\mbox{\boldmath$\bf x$},t)$ is necessary for a surrogate model with
nonlocal dynamics properties (see, e.g., [13]):
$\rho\ddot{\tilde{}\mbox{\boldmath$\bf u$}}(\mbox{\boldmath$\bf
x$})-l^{2}_{d}\nabla^{2}\tilde{}\mbox{\boldmath$\bf u$}(\mbox{\boldmath$\bf
x$})=\mbox{\boldmath$\cal L$}_{b}(\tilde{}\mbox{\boldmath$\bf
u$})(\mbox{\boldmath$\bf x$})+\mbox{\boldmath$\bf b$}(\mbox{\boldmath$\bf
x$})$, where $\rho$ is the mass density, whereas $l_{d}$ denotes an internal
dynamic scale parameter. In so doing the approach [302] is based on the
datasets ${\cal D}^{\rm DNS}=\\{\mbox{\boldmath$\bf
u$}_{k}(\mbox{\boldmath$\bf b$}_{k},\mbox{\boldmath$\bf
x$}),\mbox{\boldmath$\bf b$}_{k}(\mbox{\boldmath$\bf x$})\\}^{N}_{k=1}$
generated without compression that is the dissimilarity from the approach
(4.24)-(4.27). Some sort of compression stage is proposed in a coarse-grained,
homogenized continuum model by [303], [301] and [261] presenting a learning
framework to extract, from MD data, an optimal Linear Peridynamic Solid (LPS)
model as a surrogate for MD displacements. The coarse-graining method allows
the dependence of bond force on bond length to be determined, including the
horizon and allows substantial reductions in computational cost compared with
MD. In the original method proposed in [258], [302], the dataset ${\cal
D}^{\rm DNS}$ makes it possible to estimate the kernel $\mbox{\boldmath$\bf
K$}^{*}_{b}$ by the solution of optimization problem (4.27). However, it does
not provide reason enough to assume that the kernel $\mbox{\boldmath$\bf
K$}^{*}_{b}$ can be used for the estimation of statistical average stresses
$\mbox{$\langle$}\mbox{\boldmath$\sigma$}\mbox{$\rangle$}(\mbox{\boldmath$\bf
x$})$ (as a micromodulus was exploited for the estimation of stresses
$\mbox{\boldmath$\sigma$}(\mbox{\boldmath$\bf x$})$). Then, any information
about the statistical averages of displacements inside inclusions
$\mbox{$\langle$}\mbox{\boldmath$\bf
u$}\mbox{$\rangle$}_{i}(\mbox{\boldmath$\bf x$})$ ($\mbox{\boldmath$\bf x$}\in
v_{i}$) is lost for the online stage in the case of dataset ${\cal D}^{\rm
DNS}$, although it is precisely the operator $\mbox{\boldmath$\bf K$}^{*}_{\rm
u_{i}}$ (rather than $\mbox{\boldmath$\bf K$}^{*}_{b}$) is potentially used
for any nonlinear problem (e.g., fracture and plasticity). Moreover, we
constructed the surrogate operators (4.26) for the infinite random (or
periodic) media eliminating the edge effects and inconsistencies (see [61])
produced by popular approaches modeling the infinite inclusion fields by a
finite set of inclusions (as in [302], [304], [303], [261]). The use of body
force with compact support as a training parameter enables the edge and size
effects to be eliminated.

A new integral neural operator architecture called the Peridynamic Neural
Operator (PNO) [151] provides a surrogate operator $\cal G$
($\mbox{\boldmath$\cal G$}(\mbox{$\langle$}\mbox{\boldmath$\bf
u$}\mbox{$\rangle$})(\mbox{\boldmath$\bf x$})\approx-\mbox{\boldmath$\bf
b$}(\mbox{\boldmath$\bf x$})$) for physically consistent predictions of the
overall behavior of highly nonlinear, anisotropic, and heterogeneous materials
and achieves improved accuracy and efficiency compared to baseline models that
use predefined constitutive laws (as in Eq. (4.26)). Of course nonlineat
micromechanical model (4.15)-(4.22) can be easily incorporated into PNO by
replacement of ${\cal D}^{\rm DNS}$ (one dataset of full-field, spatial
measurements of displacement and loading) used in PNO by ${\cal D}$ (4.24) (a
few field PM datasets of compressed sets of statistical average fields)
exploited in (4.15)-(4.22). As before, it leads to the elimination of both the
edge and size effects.

However, what is even more important is that the use of ${\cal D}$ (4.24)
allows one to get a nonlocal operator counterpart of a basic concept of the LM
such as effective concentration factor (e.g.
$\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}_{i}(\mbox{\boldmath$\bf
x$})=\mbox{\boldmath$\bf A$}^{*}(\mbox{\boldmath$\bf
x$})\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}$,
$\mbox{\boldmath$\bf x$}\in v_{i}$); precisely these (“downscale”) operators
(rather than “upscale” operator $\cal G$) are potentially used for any
nonlinear problem (e.g., fracture and plasticity). It should be mentioned that
the term compressed applied to the field PM dataset ${\cal D}$ (4.24) does not
mean that the full-field dataset ${\cal D}^{\rm DNS}$ for a big sample should
be preliminary estimated. In general, the destiny of micromechanics (both the
LM and PM) is providing the required estimations by more chip, fast, robust,
and flexible methods (although it takes additional intellectual complexity to
the implementations) than DNS. On the other side, the micromechanic scheme
(4.15)-(4.22) is not a single way for estimation of the field PM datased
${\cal D}$ (4.24). So, in the computational LM, the well-known periodization
procedure [243] of random media is reduced to the modeling of periodic
boundary conditions at the unit cell containing just a few tens $n$ of
inclusions with the centers $\mbox{\boldmath$\bf x$}_{\alpha\beta}$ in
$\beta$’s Monte Carlo realization ($\alpha=1,\ldots n;\ \beta=1,\ldots,N$).
After averaging of, e.g., $\mbox{\boldmath$\sigma$}(\mbox{\boldmath$\bf
x$}_{\alpha\beta},\mbox{\boldmath$\bf x$})$ (see [35] and p. 336 in [43]) over
${\alpha,\beta}$ we can estimate
$\mbox{$\langle$}\mbox{\boldmath$\sigma$}\mbox{$\rangle$}_{i}(\mbox{\boldmath$\bf
x$})$ ($\mbox{\boldmath$\bf x$}\in v_{i}$) and ${\cal D}$ (4.24). The approach
[35] and [243] is affordable for the homogeneous (2.3) (or periodic (2.29))
boundary conditions (3.55) when the estimated
$\mbox{$\langle$}\mbox{\boldmath$\sigma$}\mbox{$\rangle$}_{i}(\mbox{\boldmath$\bf
x$})$ ($\mbox{\boldmath$\bf x$}\in v_{i}$) is translationary invariant. The
possible generalization of this approach to the inhomogeneous body force
$\mbox{\boldmath$\bf b$}(\mbox{\boldmath$\bf x$})$ requires repeating this
scheme for each $\mbox{\boldmath$\bf b$}_{k}(\mbox{\boldmath$\bf x$})$ which
is prohibitively expensive.

Comment 4.2. Effective energy-based criteria will be obtained in Section 5 for
the homogeneous VBC (2.3). Generalization of the inhomogeneous body force is
straightforward. Then, the dataset (4.24) can be extended by other effective
field parameters with a subsequent estimation of the surrogate operators for
the effective energy-based criteria.

Comment 4.3. The scheme for obtaining the set of surrogate models is presented
in Fig. 1. Block 1 DNS is used in [118], [151], [302], [304], and [303]. In
the proposed approach, we consider either the solution of Eq. (4.24) for the
random structure CMs or the fields (4.28) for each $\chi$ in the case of
periodic structure CMs. Estimations of the effective fields (4.16), (4.17),
and (4.20) (or (4.28)) are estimated in the block 2 Field PM dataset
(Compression). Solution of optimization problems (4.27) (Block 3) leads to
obtaining either one surrogate operator (e.g. $\mbox{\boldmath$\bf
K$}^{*}_{b}$ or $\cal G$, Bloc 4) or a set of surrogate models (e.g. Eq.
(4.24), Block 5). Block 2 Field PM dataset is, in fact, a union of Blocks 1-4
(without Block 5) in the Fig. 2.

Fig. 1: The scheme of obtaining of surrogate model set

###  4.4 Estimation of effective moduli

Thus, the GIE (3.32) for CM subjected to body force $\mbox{\boldmath$\bf
b$}(\mbox{\boldmath$\bf x$})$ with compact support was solved for a general
representation of the nonlinear perturbator operator $\mbox{\boldmath$\cal
L$}^{\theta\zeta}_{j}(\mbox{\boldmath$\bf
z$},\overline{\mbox{\boldmath$\zeta$}})$ (3.23). Solution (4.7) was obtained
in the framework of the closing assumption H2a) (4.1) with estimation of the
statistical averages $\mbox{$\langle$}\mbox{\boldmath$\bf
u$}\mbox{$\rangle$}_{i}(\mbox{\boldmath$\bf x$})$ (4.16) and
$\mbox{$\langle$}\mbox{\boldmath$\bf u$}\mbox{$\rangle$}(\mbox{\boldmath$\bf
x$})$ (4.19). A set of effective surrogate operators (4.26) was estimated for
the general nonlinear perturbator operator $\mbox{\boldmath$\cal
L$}^{\theta\zeta}_{q}(\mbox{\boldmath$\bf
z$},\overline{\mbox{\boldmath$\zeta$}})$ (3.23).

However, for the estimation of effective moduli (requering both the medium
statistically homogeneity and homogeneous volumetric boundary conditions
(2.31)) from the iteration approach (4.1), some additional assumptions are
needed. So, for CM with linear elastic properties of the matrix (i.e.
${}^{L}\\!\mbox{\boldmath$\bf L$}^{(0)}$ and $\mbox{\boldmath$\bf C$}^{(0)}$),
Eq. (4.22) is reduced to to the representation of the effective elastic
modulus similar to their locally elastic counterpart

$\mbox{\boldmath$\bf L$}^{*}=\mbox{\boldmath$\bf L$}^{(0)}+\mbox{\boldmath$\bf
R$}^{*},\ \ \
\mbox{$\langle$}\mbox{\boldmath$\tau$}\mbox{$\rangle$}^{l}=\mbox{\boldmath$\bf
R$}^{*}\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$},$ (4.29)

where $\mbox{\boldmath$\bf L$}^{(0)}$ is defined by Eq. (3.41) with no
assumption about $\mbox{\boldmath$\bf C$}^{(0)}$ (e.g., either the spherical
shape of ${\cal H}_{\bf x}$ or r isotropy of $\mbox{\boldmath$\bf C$}^{(0)}$).
Equation (4.292) implies linearity of the effective properties (although the
limiting representation (3.23) was obtained for generally nonlinear
perturbator operators $\mbox{\boldmath$\cal
L$}^{\theta\zeta}_{j}(\mbox{\boldmath$\bf
z$},\overline{\mbox{\boldmath$\zeta$}})$). The representation
$\mbox{$\langle$}\mbox{\boldmath$\tau$}\mbox{$\rangle$}^{l}_{i}=\mbox{$\langle$}\mbox{\boldmath$\tau$}\mbox{$\rangle$}^{l\omega}_{i}$
(4.292) can be also expressed through the displacement field at the external
interaction interface $\Gamma^{-}_{i}$ of the representative inclusion $v_{i}$
(see for details [60], [55])

$\displaystyle\mbox{$\langle$}\mbox{\boldmath$\tau$}\mbox{$\rangle$}^{l\omega}_{i}$
$\displaystyle:=$
$\displaystyle(\overline{v}_{i}^{l})^{-1}\int_{\Gamma^{-}_{i}}\mbox{\boldmath$\tau$}(\mbox{\boldmath$\bf
s$})\,d\mbox{\boldmath$\bf s$}=\mbox{$\langle$}[\mbox{\boldmath$\cal
L$}^{\sigma}(\mbox{\boldmath$\bf u$})(\mbox{\boldmath$\bf
s$})\otimes\mbox{\boldmath$\bf n$}(\mbox{\boldmath$\bf
s$})]{\,}{\,}^{\,{}^{S}}\\!\\!\\!\\!\\!\\!\otimes\mbox{\boldmath$\bf
s$}\mbox{$\rangle$}^{l\omega}_{i}$ (4.30) $\displaystyle-$
$\displaystyle\mbox{\boldmath$\bf L$}^{(0)}\mbox{$\langle$}\mbox{\boldmath$\bf
u$}(\mbox{\boldmath$\bf
s$}){\,}{\,}^{\,{}^{S}}\\!\\!\\!\\!\\!\\!\otimes\mbox{\boldmath$\bf
n$}(\mbox{\boldmath$\bf s$})\mbox{$\rangle$}^{l\omega}_{i},$

where $\mbox{\boldmath$\bf u$}(\mbox{\boldmath$\bf s$})$ ($\mbox{\boldmath$\bf
s$}\in\Gamma_{i}^{-}$) are evaluated by the use of Eqs. (3.232) and (4.7).

Subsequent simplifications of Eqs. (4.1) and (4.12) can be performed in the
framework of Hypothesis H1a (3.1) for the linear operators
$\mbox{\boldmath$\cal J$}^{\theta\zeta}_{i,j}$ (3.28) and
$\mbox{\boldmath$\cal J$}^{\theta\zeta\infty}_{i,j}$ (3.31) which can be
decomposed and reduced to the tensors at the applying to the constant
effective fields (4.1)

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\mbox{\boldmath$\cal
J$}^{I\theta\zeta}_{i,j}(\mbox{\boldmath$\bf
z$},\mbox{$\langle$}\overline{\mbox{\boldmath$\zeta$}}_{i}\mbox{$\rangle$})$
$\displaystyle+$ $\displaystyle\mbox{\boldmath$\cal
J$}^{J\theta\zeta}_{i,j}(\mbox{\boldmath$\bf
z$},\mbox{$\langle$}\overline{\mbox{\boldmath$\zeta$}}_{j}\mbox{$\rangle$})=\mbox{\boldmath$\bf
J$}^{I\theta\zeta}_{i,j}(\mbox{\boldmath$\bf
z$})\mbox{$\langle$}\overline{\mbox{\boldmath$\zeta$}}_{i}\mbox{$\rangle$}(\mbox{\boldmath$\bf
x$}_{i})+\mbox{\boldmath$\bf J$}^{J\theta\zeta}_{i,j}(\mbox{\boldmath$\bf
z$})\mbox{$\langle$}\overline{\mbox{\boldmath$\zeta$}}_{j}\mbox{$\rangle$}(\mbox{\boldmath$\bf
x$}_{j}),$ (4.31)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\mbox{\boldmath$\cal
J$}^{\theta\zeta\infty}_{i,j}(\mbox{\boldmath$\bf
z$},\mbox{$\langle$}\overline{\mbox{\boldmath$\zeta$}}_{j}\mbox{$\rangle$})$
$\displaystyle=$ $\displaystyle\mbox{\boldmath$\bf
J$}^{\theta\zeta\infty}_{i,j}(\mbox{\boldmath$\bf
z$})\mbox{$\langle$}\overline{\mbox{\boldmath$\zeta$}}_{j}\mbox{$\rangle$}(\mbox{\boldmath$\bf
x$}_{j}).$ (4.32)

In particular, for the linear operators $\mbox{\boldmath$\cal
L$}^{\theta\zeta}_{i}(\mbox{\boldmath$\bf
z$},\overline{\mbox{\boldmath$\zeta$}})$ (3.23) and homogeneous effective
strain $\overline{}\mbox{\boldmath$\varepsilon$}(\mbox{\boldmath$\bf
x$})=\overline{}\mbox{\boldmath$\varepsilon$}_{i}$ the corresponding operators
are reduced to the tensor multiplications

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\mbox{\boldmath$\vartheta$}(\mbox{\boldmath$\bf
z$})-\overline{\mbox{\boldmath$\vartheta$}}(\mbox{\boldmath$\bf z$})$
$\displaystyle=$ $\displaystyle\mbox{\boldmath$\bf
L$}^{\theta\varepsilon}_{i}(\mbox{\boldmath$\bf
z$})\overline{\mbox{\boldmath$\varepsilon$}}_{i},\ \
\mbox{\boldmath$\vartheta$}(\mbox{\boldmath$\bf z$})=\mbox{\boldmath$\bf
A$}^{\theta\varepsilon}_{i}(\mbox{\boldmath$\bf
z$})\overline{\mbox{\boldmath$\varepsilon$}}_{i},\ \
\mbox{\boldmath$\sigma$}(\mbox{\boldmath$\bf x$})=\mbox{\boldmath$\bf
L$}^{\sigma\varepsilon}_{i}(\mbox{\boldmath$\bf
x$})\overline{\mbox{\boldmath$\varepsilon$}}_{i}$ (4.33)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\mbox{\boldmath$\tau$}({\bf
x})$ $\displaystyle=$ $\displaystyle{\mbox{\boldmath$\bf
R$}}_{i}(\mbox{\boldmath$\bf x$})\overline{\mbox{\boldmath$\varepsilon$}}_{,}$
(4.34)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\mbox{\boldmath$\bf
R$}_{i}(\mbox{\boldmath$\bf
x$})\mbox{$\langle$}\overline{\mbox{\boldmath$\varepsilon$}}\mbox{$\rangle$}$
$\displaystyle=$ $\displaystyle\mbox{\boldmath$\cal
L$}^{\sigma}(\mbox{\boldmath$\bf
C$}_{1}\mbox{\boldmath$\eta$})(\mbox{\boldmath$\bf x$}),\ \
\mbox{\boldmath$\eta$}(\mbox{\boldmath$\bf z$},\mbox{\boldmath$\bf
y$})=\big{[}\mbox{\boldmath$\bf A$}^{u\epsilon}_{i}(\mbox{\boldmath$\bf
z$}-\mbox{\boldmath$\bf x$}_{i})-\mbox{\boldmath$\bf
A$}^{u\epsilon}_{i}(\mbox{\boldmath$\bf y$}-\mbox{\boldmath$\bf
x$}_{i})\big{]}\overline{\mbox{\boldmath$\varepsilon$}}.$ (4.35)

Then, substitutions of Eqs. (4.31) and (4.32) into Eq. (4.1) leads to the
linear algebraic equation

$\displaystyle\langle\overline{\mbox{\boldmath$\vartheta$}}_{i}\rangle(\mbox{\boldmath$\bf
z$})$ $\displaystyle=$
$\displaystyle\langle\mbox{\boldmath$\vartheta$}\rangle({\bf
z})+\int\big{\\{}[\mbox{\boldmath$\bf
J$}^{I\theta\zeta}_{i,j}(\mbox{\boldmath$\bf
z$})\mbox{$\langle$}\overline{\mbox{\boldmath$\zeta$}}_{i}\mbox{$\rangle$}(\mbox{\boldmath$\bf
x$}_{i})+\mbox{\boldmath$\bf J$}^{J\theta\zeta}_{i,j}(\mbox{\boldmath$\bf
z$})\mbox{$\langle$}\overline{\mbox{\boldmath$\zeta$}}_{j}\mbox{$\rangle$}(\mbox{\boldmath$\bf
x$}_{j})]$ (4.36) $\displaystyle\times$ $\displaystyle\varphi(v_{j},{\bf
x}_{j}|;v_{i},{\bf x}_{i})-\mbox{\boldmath$\bf
J$}_{i,j}^{\theta\zeta\infty}(\mbox{\boldmath$\bf
z$})\mbox{$\langle$}\overline{\mbox{\boldmath$\zeta$}}_{j}\mbox{$\rangle$}(\mbox{\boldmath$\bf
x$}_{j})\bigl{\\}}d{\bf x}_{j}.$

Equations (4.36) and (4.34) can be explicitly solved with the estimation of
the effective tensors $\mbox{\boldmath$\bf R$}^{*}$ (4.292) and the effective
moduli $\mbox{\boldmath$\bf L$}^{*}$ (4.29). This scheme was realized for the
linear bond-based model of constituents of CM [57], and for the linear state-
based thermoelastic model [64], [67]. In particular, the effective field
method (EFM) used the closing assumptions H2b (4.2) and (4.3)

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\mbox{\boldmath$\bf
L$}^{\rm*EFM}$ $\displaystyle=$ $\displaystyle\mbox{\boldmath$\bf
L$}^{(0)}+c^{l(1)}\mbox{\boldmath$\bf R$}_{i}^{l\omega}\mbox{\boldmath$\bf
D$}^{\rm EFM},$ (4.37)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!(\mbox{\boldmath$\bf
D$}^{\rm EFM})^{-1}$ $\displaystyle=$ $\displaystyle\mbox{\boldmath$\bf
I$}-\int\mbox{$\langle$}\mbox{\boldmath$\bf n$}(\mbox{\boldmath$\bf
s$}){\,}{\,}^{\,{}^{S}}\\!\\!\\!\\!\\!\\!\otimes\mbox{\boldmath$\bf
L$}_{j}^{u\epsilon}(\mbox{\boldmath$\bf s$}-\mbox{\boldmath$\bf
x$}_{i})\mbox{$\rangle$}^{l\omega}_{i}[\varphi(v_{j},{\bf x}_{j}|;v_{i},{\bf
x}_{i})-n^{(j)}]d{\bf x}_{j},$ (4.38)

where $\mbox{\boldmath$\bf R$}_{i}^{l\omega}=\mbox{\boldmath$\bf R$}_{i}^{l}$.
For the Mori-Tanaka approach (MTA) using the closing assumption

$\displaystyle\overline{\mbox{\boldmath$\varepsilon$}}(\mbox{\boldmath$\bf
x$})=\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}^{(0)}$
(4.39)

instead of the assumptions (4.2) and (4.3), we get $(\mbox{\boldmath$\bf
s$}\in v_{i}^{l}$)

$\displaystyle\mbox{\boldmath$\bf L$}^{\rm*MTA}$ $\displaystyle=$
$\displaystyle\mbox{\boldmath$\bf L$}^{(0)}+c^{l(1)}\mbox{\boldmath$\bf
R$}_{i}^{l\omega}\mbox{\boldmath$\bf D$}^{\rm MTA},$ (4.40)
$\displaystyle(\mbox{\boldmath$\bf D$}^{\rm MTA})^{-1}$ $\displaystyle=$
$\displaystyle c^{l(0)}\mbox{\boldmath$\bf I$}+c^{l(1)}\mbox{\boldmath$\bf
A$}_{i}^{l\omega},\ \ \mbox{\boldmath$\bf
A$}_{i}^{l\omega}:=\mbox{$\langle$}\mbox{\boldmath$\bf
A$}^{u\epsilon}_{i}(\mbox{\boldmath$\bf s$})\mbox{$\rangle$}^{l\omega}_{i}.$
(4.41)

Numerical results for 1D cases are presented in [47], [48], [50], [57].

Comment 4.4. For local elasticity, this method (called the second background
of micromechanics) uses the new GIE proposed in [44], [45] (see for details
[60]) and permits to abandonment of the basic concepts of micromechanics: the
hypothesis of “ellipsoidal symmetry”, and the effective field hypothesis EFH,
H1a (3.1). The results of this abandonment for local elasticity was
quantitatively estimated for some models of CMs reinforced by aligned
continuously inhomogeneous fibers. Some new effects (with possible correction
of a sign of statistical averages of local field estimations) were discovered
that were impossible in the framework of the classical (the first) background
of micromechanics. Representations (4.23)-(4.26) (with the corresponding
estimations of effective material and field parameters) can be considered as
generalizations of the mentioned linear models of LM to their peridynamic
counterparts. Furthermore, the LM methods mentioned are linear whereas the PM
methods (see Eqs. (3.23)-(3.36), (4.1)-(4.1), and (4.22)) are nonlinear. We
used only linearity for the matrix at the obtaining of Eqs. (4.17), (4.171),
and (4.18). It would be interesting to obtain the nonlinear LM counterparts of
the mentioned method in PM. I.e. we have a rare (the second) case when the
methods developed in PM initiate the new methods of LM unknown before.

Comment 4.5. It should be mentioned that the PM counterpart of the LM method’s
classification (3.1). So, consistent4.5)-(4.19) corresponts to the LM self-
consisten methods Gr4) (3.1). Perturbation method Gr2) corresponds to the
dilute approximation in the PM considered in [58]. The popular method of the
effective medium (MEM) belonging to class Gr1) in the LM has a straightforward
generalization to the PM as one peridynamic inclusion in the local effective
medium. Variational methods were considered in [59]: the theorem of work and
energy, Betti’s reciprocal theorem, and the theorem of virtual work were
proved. Principles of a minimum of both potential energy and complementary
energy were generalized. It was obtained the generalized Hill’s bounds on the
effective elastic moduli of peridynamic random structure composites. Of
special interest in the LM are Hashin-Shtrilman (H-S) bounds, which were, in
particular, obtained from the GIE (3.4), see [289]. However, H-S bounds were
also estimated from the GIE (3.6) (see pp. 439-442 in [60]) or CMs with
noncanonical inhomogeneous inclusions; the difference between the H-S bounds
obtained from Eqs. (3.4) and (3.6) is described by a 2.9. Owing to the formal
similarity of GIEs in the LM (3.6), (3.10) and PM (3.34), it would be
interesting to generalize the H-S bounds obtained from Eq. (3.6) to their
peridynamic counterpart obtained from the GIE (3.34).

###  4.5 Damage mechanics and microcracked media

Of special interest is damage mechanics. Damage is a basic concept in
engineering design exploited for predictions about the initiation,
propagation, and fracture of materials without molecular scale, which can be
too complex for engineering analysis (see the pioneering works by Kachanov
[155] and Rabotnov [237] as well as [156], [168]). Physically justified damage
theory was proposed in peridynamics (see [254], [259]), where damage is
introduced by permitting damaged bonds to break irreversibly, i.e., damaging
the bond is considered as removing if the bond stretch exceeds a critical
stretch for some time. In other words, the damaged medium is considered as a
composite material in which some two-point micromodulus functions vanish. We
estimate the effective moduli of the damaged medium (initially homogeneous) as
the functions of damage, damage functions, and micromodulus profile. In so
doing, all mentioned parameters of damage are fixed and assumed to be known,
i.e., we don’t consider the accumulation of damage and their dependence on the
applied loading as well as crack initiation and propagation (see, e.g., [3],
[134]).

Following [259], we introduce a damage function $\mu(\mbox{\boldmath$\bf
x$},\mbox{\boldmath$\xi$})$ and local damage $\varphi(\mbox{\boldmath$\bf
x$})$at a point $\bf x$

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\mu(\mbox{\boldmath$\bf
x$},\mbox{\boldmath$\xi$})=\begin{cases}0,\\!&{\rm the\ bond}\
\mbox{\boldmath$\xi$}\ {\rm is\ missing},\\\ 1,\\!&{\rm the\ bond}\
\mbox{\boldmath$\xi$}\ {\rm is\ acting},\end{cases}\\!\\!\\!,\ \
\varphi(\mbox{\boldmath$\bf x$})=1-\\!\\!\int_{{\cal H}_{\bf
x}}\\!\\!\\!\mu(\mbox{\boldmath$\bf
x$},\mbox{\boldmath$\xi$})d\mbox{\boldmath$\xi$}/\\!\\!\int_{{\cal H}_{\bf
x}}\\!\\!\\!d\mbox{\boldmath$\xi$}.$ (4.42)

Then for an undamaged (pristine) material $\mbox{\boldmath$\bf
C$}(\mbox{\boldmath$\bf x$},\hat{}\mbox{\boldmath$\bf
x$})\equiv\mbox{\boldmath$\bf C$}^{(0)}(\mbox{\boldmath$\bf
x$},\hat{}\mbox{\boldmath$\bf x$})$ while the micromodulus of damaged medium
can be presented in the form $\mbox{\boldmath$\bf C$}(\mbox{\boldmath$\bf
x$},\hat{}\mbox{\boldmath$\bf x$})=\mu(\mbox{\boldmath$\bf
x$},\hat{}\mbox{\boldmath$\bf x$})\mbox{\boldmath$\bf
C$}^{(0)}(\mbox{\boldmath$\bf x$},\hat{}\mbox{\boldmath$\bf x$})$. We
decompose this micromodulus as $\mbox{\boldmath$\bf C$}(\mbox{\boldmath$\bf
x$},\hat{}\mbox{\boldmath$\bf x$})=\mbox{\boldmath$\bf
C$}^{(0)}(\mbox{\boldmath$\bf x$},\hat{}\mbox{\boldmath$\bf
x$})+\mbox{\boldmath$\bf C$}_{1}(\mbox{\boldmath$\bf
x$},\hat{}\mbox{\boldmath$\bf x$})$, where $\mbox{\boldmath$\bf
C$}_{1}(\mbox{\boldmath$\bf x$},\hat{}\mbox{\boldmath$\bf
x$})=(\mu(\mbox{\boldmath$\bf x$},\mbox{\boldmath$\xi$})-1)\mbox{\boldmath$\bf
C$}^{(0)}(\mbox{\boldmath$\bf x$},\hat{}\mbox{\boldmath$\bf x$})$ is a jump of
micromodulus with respect to the micromodulus of the matrix (undamaged
medium). Thus, we reformulated the damage model in the term of PM, which was
solved, e.g., for periodic media (see for details [60]), when we estimate
$\mbox{$\langle$}\mbox{\boldmath$\sigma$}\mbox{$\rangle$}^{\Omega}$,
$\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}^{\Omega}$, and
$\mbox{\boldmath$\bf L$}^{*}$. It was proved that for uniform damage functions
$\mu(\mbox{\boldmath$\bf x$},\mbox{\boldmath$\xi$})$, for any profile either
(2.251) and (2.252), and any dimension $d=1,2,3$, $\mbox{\boldmath$\bf
L$}^{*}$ is close to the popular engineering approximation (see, e.g. [156])
$\mbox{\boldmath$\bf L$}^{*}=[1-\varphi(\mbox{\boldmath$\bf
x$})]\mbox{\boldmath$\bf L$}^{(0)}$.

It should be mentioned that the concept of peridynamic damage (initiated by
Silling [254], see also [259]) expressed through the damage function
$\mu(\mbox{\boldmath$\bf x$},\mbox{\boldmath$\xi$})$ is universally recognized
by the Peridynamic Society and widely used for analyses of the finite-size
samples of either initially homogeneous media or CMs. However, the infinite
damaged media is of special interest just because of the classical continuum
damage (CD) theory (see the books [156], [168]) was at first conceptually
introduced for statistically homogeneous media as a state variable affecting
the exploitation properties of this infinite media. In other words, so-called
the scale separation hypothesis is implicitly assumed to be valid, which means
that both the sample scale, and the applied field scale infinitely exceed a
damage scale. Moreover, it is also interesting to set up relationships between
the CD theories and peridynamic damage.

The popular idea for the incorporation of existing classical CD models into
the peridynamic damage theory is an approximation of some peridynamic
parameters by a function depending on a single parameter $\varphi$, which is
commonly the case in some of the most popular damage models used in practice.
So, one [277] proposed to describe an influence function of damaged medium to
be proportional to the conventional (radial) influence function (for a
pristine medium) with the proportionality factor $(1-\varphi)$ depending on
the state of damage between two points defining a bond. Han et al. [135] used
an approximation for the damaged stiffness tensor (called the effective
stiffness in our notation) defined by the micromodulus, which in turn depends
on the critical local damage rather than on the damage function. Thus, the
known incorporations of the CD concept into peridynamics are performed by the
presentations of one or another peridynamic parameters (e.g. either influence
function, micropotential, or micromodulus) as the functions of a state of
damage (which is commonly the case in some of the most popular damage models
used in practice). In the fundamental work by Littlewood et al. [181], the
damage state (which is conceptually more general than the local damage
$\varphi$) is introduced in the most general thermodynamic framework with the
strain energy density function presented as a dot product of micropotential of
undamaged bonds, and the damage state.

Numerical results were obtained for 1D case for the different damage function
$\mu(\mbox{\boldmath$\bf x$},\mbox{\boldmath$\xi$})$ and the different
micromodulus profiles (2.251) and (2.252) (see for details [54], [60]). It was
demonstrated that for fixed local damage $\varphi=0.5$, the effective modulus
$E^{*}$ has significantly more sensitive dependence on details of both the
damage parameters (such as the damage function $\mu(\mbox{\boldmath$\bf
x$},\mbox{\boldmath$\xi$})$) and the peristatic properties (such as
micromodulus profile either (2.251) and (2.252)) rather than on only one
parameter of the local damage $\varphi$. Only in a limiting case of both the
homogeneous distribution of damaged bonds within the horizon, and the large
node number $n^{l}$, the peridynamic periodic damage model is reduced to the
classical CD model $\mbox{\boldmath$\bf L$}^{*}=[1-\varphi(\mbox{\boldmath$\bf
x$})]\mbox{\boldmath$\bf L$}^{(0)}$ [156]. In some sense, a goal of the
proposed approach can be considered as an inversion with respect to the goals
of the known approaches (see [135], [277], and the most general version of
peridynamic damage theory proposed in [181]). We don’t try to incorporate the
concept of the CD theory into peridynamics. Rather the opposite, we extracted,
in particular, apparently the most popular equation of CD theory
$\mbox{\boldmath$\bf L$}^{*}=[1-\varphi(\mbox{\boldmath$\bf
x$})]\mbox{\boldmath$\bf L$}^{(0)}$ as a particular case of the most simple
peridynamic damage theory (4.27). Moreover, in the framework of the functional
graded theory (see for details [60]) it was proposed the model for the
peristatic generalization of the classical feature of damage localization. The
detected strain redistributions in damaged material quantitatively explain a
qualitative feature of damage proposed in [12] when the breakage of individual
bonds more likely generates its neighbors break, because more load is shifted
to their neighbors.

In local micromechanics, the determination of effective stiffness for bodies
containing microcracks is one of the fundamental aspects of micromechanical
modeling. This subject is usually treated from various theoretical standpoints
corresponding to one or another class of methods developed for micromechanics
of CMs (3.1) (see the reviews, e.g., [43], [60], [223], [275], [289]). Among
the self-consistent group, apparently the effective medium method [169],
[140], Mori-Tanaka approach [216], and effective field method [287] (see also
[160], [162]) are most popular and used for their subsequent generalizations
(see for references, e.g. [103], [60]). So, the self-consistent method (some
sort of effective medium method) by [25] consists of the calculation of crack
opening displacements (COD) for microcracks embedded in the unknown effective
medium. When an infinite number of randomly distributed inhomogeneities are
analyzed, the average field in the matrix can be considered as a good
approximation of the effective field that was pointed out by Benveniste [26]
in a reformulation of the MTM which was also applied to the cracked medium
[22]; its popularity is probably based upon the simplicity of its application.
The effective field hypothesis (EFH, also called the H1a hypothesis, p. 253 in
[43], see Section 3) is the most fundamental and most exploited concept of
micromechanics. It has been demonstrated [284] that the MT approach reduces to
the formula by Willis [287] (see also [236]) in the case of a matrix
containing a set of aligned ellipsoidal inclusions, and this, apparently, is
considered a virtue. For the flat distribution of cracks, the formulae by
Willis [287] and MT approach [22] are identical because the MT procedure
implicitly assumes the same distribution for the inclusions as their shape,
which corresponds to the general estimate [236] for the case where the
distribution and inclusion have the same aspect ratio.

Of course, a microcracked medium can be considered as a composite material
where $\mbox{\boldmath$\bf C$}^{(1)}(\mbox{\boldmath$\bf
x$},\hat{}\mbox{\boldmath$\bf x$})=-\mbox{\boldmath$\bf
C$}^{(0)}(\mbox{\boldmath$\bf x$},\hat{}\mbox{\boldmath$\bf x$})$ (as
${}^{L}\\!\mbox{\boldmath$\bf L$}^{(1)}(\mbox{\boldmath$\bf
x$})=-^{L}\\!\mbox{\boldmath$\bf L$}^{(0)}$ and all methods of the PM
considered in Subsection 4.1-4.4 are applicable for such medium.

###  4.6 Schematic representation of nonlinear GIE implementation

The scheme of the CAM (3.23), (3.28), (4.1)-(4.8) is presented in the blok
structure in Fig. 2. For shortening of notations, the presented scheme
corresponds to the iteration solution of

Fig. 2: The iteration scheme of the CAM (4.19)

Eq. (4.1). The central ingredient contains the Blocks 3 Perturbator and 4
Micromechanics placed in a circle. The perturbator $\mbox{\boldmath$\cal
L$}^{\theta\zeta}_{i}(\mbox{\boldmath$\bf
z$},\overline{\mbox{\boldmath$\zeta$}})$ (3.28) (and related perturbators is
evaluated in Block 3. A preprocessing 1 Input to Block 3 contains geometrical
microtopology of one inclusion $v_{i}(\mbox{\boldmath$\bf x$})$ in the
infinite homogeneous matrix (or large sample), material properties
$\mbox{$\mathbb{C}$}(\mbox{\boldmath$\bf x$})$ (e.g., micromodulus
$\mbox{\boldmath$\bf C$}^{\rm bond}$ (2.21)), and square mesh $\Omega^{\rm
sq}:=\\{(x_{1},...,x_{d})^{\top}|x_{i}=hp_{i}\\}$ $(\mbox{\boldmath$\bf
p$}=(p_{1},\ldots,p_{d})\in Z^{d})$. A preprocessing 2 Input to Block 4 only
contains geometrical information: the same mesh $\Omega^{\rm sq}$ and the
distribution of inclusions $X^{\partial v}$ (described by the probability
densities $\varphi(v_{j},{\bf x}_{j})$ and $\varphi(v_{j},{\bf
x}_{j}|;v_{i},{\bf x}_{i})$ for the random structures). The need to use the
same mesh $\Omega^{\rm sq}$ in Block 3 and Block 4, is explained by the
convenience of matching these Blocks (see the case of LM on p. 428 in [60]).
In the cycle [Block 3]$\rightleftharpoons$[Block 4] (4.1), the input-output of
the Block 3 are
$\mbox{$\langle$}\overline{\mbox{\boldmath$\zeta$}}\mbox{$\rangle$}_{i}(\mbox{\boldmath$\bf
x$})$ and $\mbox{\boldmath$\cal L$}^{\theta\zeta}_{i}(\mbox{\boldmath$\bf
z$},\mbox{$\langle$}\overline{\mbox{\boldmath$\zeta$}}\mbox{$\rangle$}_{i})$ ,
respectively, where $\mbox{\boldmath$\bf x$},\mbox{\boldmath$\bf
z$}\in\Omega^{\rm sq}$. I.e. the mesh $\Omega^{\rm sq}$ is only used for
compatibility of Blocks 3 and 4 rather than for a solution of the problem
inside Block 3 (see Eqs. (3.23)). Any inside knowledge of Bloc 3 (or 4) is not
required in Blok 4 (or 3), respectively. Therefore, Blocks 3 and 4 can be
considered as the Black Boxes (visualized in Fig. 2 by black color) to each
other.

In Fig. 3, a zoomed Block 3 Perturbator is presented for one step in the
circle [Block 3]$\rightleftharpoons$[Block 4] (see Fig.2). The Input Block in
Fig. 3 presents the effective field
$\mbox{$\langle$}\overline{\mbox{\boldmath$\vartheta$}}\mbox{$\rangle$}_{i}(\mbox{\boldmath$\bf
x$})$ (estimated at the previous step in Block 4, see Eq. (4.22)) acting in
the vicinity of the inclusion $v_{i}$ (4.3)) replaced by the matrix material.
The Output block in Fig. 3 contains the perturbation field
$\mbox{$\langle$}\mbox{\boldmath$\vartheta$}\mbox{$\rangle$}_{i}(\mbox{\boldmath$\bf
z$})-\mbox{$\langle$}\overline{\mbox{\boldmath$\vartheta$}}\mbox{$\rangle$}_{i}(\mbox{\boldmath$\bf
z$})=\mbox{\boldmath$\cal L$}^{\theta\zeta}_{i}(\mbox{\boldmath$\bf
z$},\mbox{$\langle$}\overline{\mbox{\boldmath$\zeta$}}\mbox{$\rangle$}_{i})$
(4.3) ($\mbox{\boldmath$\bf z$}\in\mathbb{R}^{d}$) inside and outside the
inclusion $v_{i}$ in the local square mesh $\Omega^{\rm sq}_{\rm local}$
assigned to the inclusions $v_{i}$. In so doing, the mesh used for the
solution of Eq. (4.3) (e.g., FEA mesh, an internal mesh of Block 3, which is
not presented in Fig. 3) has no connection with the mesh $\Omega^{\rm sq}_{\rm
local}$ exploited as an external mesh to Block 3 and utilized only for the
interaction of Block 3 and Block 4.

Fig. 3. Zoomed Block 3 (from Fig. 2.) with the input effective field and
output perturbation field

Zoomed Block 4 is presented in Fig. 4, where the square mesh $\Omega^{\rm sq}$
has two

Fig. 4: Zoomed Block 4 (from Fig. 2.) with the input perturbation field and
output effective field

constitutive meshes. The moving local mesh $\Omega^{\rm sq}_{\rm local}$
(solid black grid) is assigned to the moving inclusion $v_{j}$ (see Eq.
(4.19)) with the center $\mbox{\boldmath$\bf x$}_{j}$ so that
$\mbox{\boldmath$\bf x$}_{j}={\bf 0}$ in the system $\Omega^{\rm sq}_{\rm
local}$. An effect zone (covered by the local mesh $\Omega^{\rm sq}_{\rm
local}$) of field perturbations produced by the heterogeneity $v_{j}^{l}$ is
limited by the Minkowski addition of domains $\mbox{\boldmath$\bf x$}\in
v_{j}^{\oplus}\equiv v_{j}^{\phi}\oplus v_{j}^{l}$ (see for details [57]),
where $v_{j}^{\phi}$ is a domain for location of inclusion centers
$\mbox{\boldmath$\bf x$}_{j}$ were $\varphi(v_{j},\mbox{\boldmath$\bf
x$}_{j}|v_{i},\mbox{\boldmath$\bf x$}_{i})-n^{(j)}$ is not negligible (in the
sense of influence on $\mbox{\boldmath$\bf L$}^{*}$). The global mesh
$\Omega^{\rm sq}_{\rm global}$ (dashed black grid) is assigned to the fixed
inclusion $v_{i}$ (see Eq. (4.19)) with the center $\mbox{\boldmath$\bf
x$}_{i}$ so that $\mbox{\boldmath$\bf x$}_{i}={\bf 0}$ in the system
$\Omega^{\rm sq}_{\rm global}$ and $\mbox{\boldmath$\bf x$}_{j}\in\Omega^{\rm
sq}_{\rm global}$ ($j=1,2,\ldots$). The main reason for the choice of the same
square mesh $\Omega^{\rm sq}_{\rm local}$ and $\Omega^{\rm sq}_{\rm global}$
is explained in [57]. For statistically homogeneous CMs subjected to the
homogeneous VBC (2.3), the input perturbator $\mbox{\boldmath$\cal
L$}^{\theta\zeta}_{j}(\mbox{\boldmath$\bf
z$},\mbox{$\langle$}\overline{\mbox{\boldmath$\zeta$}}\mbox{$\rangle$}_{j})$
($\mbox{\boldmath$\bf z$}\in\mathbb{R}^{d}$) assigned to the moving inclusion
$v_{j}$ ($j=1,2,\ldots$) is estimated only one time. The output of the Block 4
is the effective field
$\mbox{$\langle$}\overline{\mbox{\boldmath$\zeta$}}\mbox{$\rangle$}_{i}(\mbox{\boldmath$\bf
x$})$ acting on the fixed inclusion $v_{i}$.

Zoomed Block 5 from Fig. 2 is presented in Fig. 5. Here the effective field
$\mbox{$\langle$}\overline{\mbox{\boldmath$\zeta$}}\mbox{$\rangle$}_{i}(\mbox{\boldmath$\bf
x$})$ ($\mbox{\boldmath$\bf x$}\in v_{i}$) acting on the representative area
of the inclusion $v_{i}$ (replaced by the matrix material) are estimated as
output of the Block 4 (see Fig. 2) at the last step in the circle [Block
3]$\rightleftharpoons$[Block 4]. The field inside the representative inclusion
$v_{i}$
$\mbox{$\langle$}\mbox{\boldmath$\vartheta$}\mbox{$\rangle$}_{i}(\mbox{\boldmath$\bf
z$})=\mbox{$\langle$}\overline{\mbox{\boldmath$\vartheta$}}\mbox{$\rangle$}_{j}(\mbox{\boldmath$\bf
z$})+\mbox{\boldmath$\cal L$}^{\theta\zeta}_{i}(\mbox{\boldmath$\bf
z$},\overline{\mbox{\boldmath$\zeta$}})$ ($\mbox{\boldmath$\bf z$}\in v_{i}$)
is estimated by the use of the Block 3 (see Fig. 2) at the last step in the
circle [Block 3]$\rightleftharpoons$[Block 4]. At last, the effective
constitutive equation (4.23) is evaluated only one time through the
substitution of
$\mbox{$\langle$}\mbox{\boldmath$\vartheta$}\mbox{$\rangle$}_{j}(\mbox{\boldmath$\bf
z$})$ into Eqs. (4.21). In the case of the homogeneous VBC (2.3), the
effective moduli $\mbox{\boldmath$\bf L$}^{*}$ (4.29) are expressed (also only
one time) through the statistical average displacement field
$\mbox{$\langle$}\mbox{\boldmath$\bf u$}(\mbox{\boldmath$\bf
u$})\mbox{$\rangle$}_{i}(\mbox{\boldmath$\bf s$})$ at the external interaction
interface ($\mbox{\boldmath$\bf s$}\in\Gamma_{i}$ marked by a black color) of
the representative heterogeneity $v_{i}$; i.e. the field
$\mbox{$\langle$}\mbox{\boldmath$\vartheta$}\mbox{$\rangle$}_{i}(\mbox{\boldmath$\bf
z$})$ incide the inclusion $\mbox{\boldmath$\bf z$}\in v_{i}$ is bot used at
the estimation for estimation of $\mbox{\boldmath$\bf L$}^{*}$.

Fig. 5: Zoomed Block 5 (from Fig. 2.) with the field
$\mbox{$\langle$}\mbox{\boldmath$\vartheta$}\mbox{$\rangle$}_{i}(\mbox{\boldmath$\bf
z$})$ (dashed line) inside inclusion $v_{i}$

###  4.7 Nonlinear incremental Mori-Tanaka approach in PM

In general, within a FEA ABAQUS specifies strain increments, which are taken
as the input for the constitutive material law. The stress response and the
consistent (algorithmic) tangent tensor of the composite material at the end
of the increment have to be evaluated for this strain increment. The following
strategy of the incremental Mori-Tanaka approach is used for solving the
nonlinear behavior of peridynamic CM (which is a generalization of
corresponding approaches in the LM, see, e.g., [100], [108], [161], [233]).

At first, we consider a macroscopic matrix domain $w$ subjected to the remote
homogeneous volumetric boundary conditions (2.3) producing a linear
displacement field $\mbox{\boldmath$\bf u$}(\mbox{\boldmath$\bf
x$})=\mbox{\boldmath$\varepsilon$}^{w\Gamma}\cdot\mbox{\boldmath$\bf x$}$
($\mbox{\boldmath$\bf x$}\in w)$. For nonlinear PD, the tangent moduli
$\mbox{\boldmath$\bf L$}^{(0)}_{t}$ are represented and approximated by a
forward difference of the associated stresses

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\Delta\mbox{\boldmath$\sigma$}$
$\displaystyle=$ $\displaystyle\mbox{\boldmath$\bf
L$}^{(0)}_{t}\Delta\mbox{\boldmath$\varepsilon$}^{w\Gamma},$ (4.43)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\mbox{\boldmath$\bf
L$}^{(0)}_{t}\Delta\mbox{\boldmath$\varepsilon$}^{w\Gamma}$ $\displaystyle=$
$\displaystyle\big{[}\mbox{\boldmath$\cal
L$}^{\sigma(0)}((\mbox{\boldmath$\varepsilon$}^{w\Gamma}+\Delta\mbox{\boldmath$\varepsilon$}^{w\Gamma})\cdot\mbox{\boldmath$\bf
x$})-\mbox{\boldmath$\cal
L$}^{\sigma(0)}(\mbox{\boldmath$\varepsilon$}^{w\Gamma}\cdot\mbox{\boldmath$\bf
x$})\big{]}\Delta\mbox{\boldmath$\varepsilon$}^{w\Gamma},$ (4.44)

where $\mbox{\boldmath$\bf L$}^{(0)}_{t}$ represents the associated tangent
moduli, or the components of the tangent modulus tensor obtained by the
perturbation of $\Delta\mbox{\boldmath$\varepsilon$}^{w\Gamma}$.

We assume that up to $n$th increment
$\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}^{[n]}$ the
behavior of CM is linear and we can estimate $\mbox{\boldmath$\bf
L$}^{(0)[n]}$ (3.44), $\mbox{\boldmath$\bf L$}^{u\epsilon[n]}$ (4.331),
$\mbox{\boldmath$\bf L$}^{\sigma\epsilon[n]}$ (4.342), $\mbox{\boldmath$\bf
R$}_{i}^{l\omega[n]}$ (4.341), $\mbox{\boldmath$\bf A$}_{i}^{l\omega[n]}$
(4.332),
$\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}^{l(0)[n]}\equiv\overline{}\mbox{\boldmath$\varepsilon$}^{l(0)[n]}$
(4.39), $\mbox{\boldmath$\bf L$}_{i}^{*\omega{\rm MTA}[n]}$ (4.40). In a
backward Euler fully implicit numerical scheme, the value of a given quantity
$\chi$ is updated from the previous time step $n$ to the current $n+1$ per
$\mbox{\boldmath$\chi$}^{[n+1]}=\mbox{\boldmath$\chi$}^{[n]}+\Delta\mbox{\boldmath$\chi$}^{[n+1]}$.
Such an implicit relation is usually solved iteratively, and the current value
is updated for each iteration by:
$\mbox{\boldmath$\chi$}^{[n+1](k+1)}=\mbox{\boldmath$\chi$}^{[n+1](k)}+\Delta\mbox{\boldmath$\chi$}^{[n+1](k)}$
until $\mbox{\boldmath$\chi$}^{[n+1](k+1)}$ has converged. x(n+1)(k).

Using the Backward Euler framework, equation (4.27) is written in incremental
form at time step $n+1$ as

$\displaystyle\Delta\mbox{$\langle$}\mbox{\boldmath$\sigma$}\mbox{$\rangle$}^{[n+1]}=\mbox{\boldmath$\bf
L$}^{*\omega{\rm
MTM}[n+1]}\Delta\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}^{[n+1]},$
(4.45)

The known global state and instantaneous behavior from the previous increment
together with the prescribed strain increment,
$\Delta\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}^{[n+1]}$,
are used to calculate a prediction for the new ‘homogenized’ overall state.
This predicted new state gives rise to a predicted matrix strain increment
$\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}^{l(0)[n]}$,
which is used to solve for matrix perturbators $\mbox{\boldmath$\bf
L$}^{u\epsilon[n]}$ (4.331) and $\mbox{\boldmath$\bf L$}^{\sigma\epsilon[n]}$
(4.342) in an Euler backward scheme on the micro level. In an incremental-
iterative form, the prediction step states that the stress during the current
loading increment $n+1$ and at the end of each iteration $k+1$ is given as

$\displaystyle\mbox{$\langle$}\mbox{\boldmath$\sigma$}\mbox{$\rangle$}^{[n+1](k+1)}=\mbox{$\langle$}\mbox{\boldmath$\sigma$}\mbox{$\rangle$}^{[n]}+\mbox{\boldmath$\bf
L$}^{*\omega{\rm
MTM}[n+1](k+1)}\Delta\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}^{[n+1]},$
(4.46)

where $\mbox{\boldmath$\bf L$}^{*\omega{\rm MTM}[n+1](1)}=\mbox{\boldmath$\bf
L$}^{*\omega{\rm MTM}[n]}$ and the initial iteration corresponds to $k=1$,
i.e. $\mbox{\boldmath$\bf L$}^{(0)[n](1)}=\mbox{\boldmath$\bf L$}^{(0)[n]}$
(3.44), $\mbox{\boldmath$\bf L$}^{u\epsilon[n](1)}=\mbox{\boldmath$\bf
L$}^{u\epsilon[n]}$ (4.331), $\mbox{\boldmath$\bf
L$}^{\sigma\epsilon[n](1)}=\mbox{\boldmath$\bf L$}^{\sigma\epsilon[n]}$
(4.342)),
$\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}^{l(0)[n](1)}=\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}^{l(0)[n]}$,
$\mbox{\boldmath$\bf R$}_{i}^{l\omega[n](1)}=\mbox{\boldmath$\bf
R$}_{i}^{l\omega[n]}$ (ref4.341), $\mbox{\boldmath$\bf
A$}_{i}^{l\omega[n](1)}=\mbox{\boldmath$\bf A$}_{i}^{l\omega[n]}$ (4.332),
$\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}^{l(0)[n](1)}=\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}^{l(0)[n]}$
(4.39), $\mbox{\boldmath$\bf L$}_{i}^{*\omega{\rm
MTM}[n](1)}=\mbox{\boldmath$\bf L$}_{i}^{*\omega{\rm MTM}[n]}$ (4.40).

For the current
$\Delta\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}^{[n+1]}$
we estimate
$\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}^{l(0)[n+1](k)}$
by Eq. (4.40) with $\mbox{\boldmath$\bf A$}_{i}^{l\omega}=\mbox{\boldmath$\bf
A$}_{i}^{l\omega(k)}$.
$\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}^{l(0)[n+1](k)}$
is substituted to the right-hand side of Eqs. (4.43) , (4.331) and (4.342)
with estimations of $\mbox{\boldmath$\bf L$}_{t}^{(0)(k+1)}$ (4.44),
$\mbox{\boldmath$\bf L$}^{u\epsilon[n](k+1)}$ (4.331), and
$\mbox{\boldmath$\bf L$}^{\sigma\epsilon[n](k+1)}$ (4.342) which, in its turn,
prove estimations of $\mbox{\boldmath$\bf R$}_{i}^{l\omega[n](k+1)}$ (4.40),
$\mbox{\boldmath$\bf A$}_{i}^{l\omega[n](k+1}$ (4.41), and
$\mbox{\boldmath$\bf L$}_{i}^{*\omega{\rm MTM}[n](k+1)}$ (4.40). The process
is completed when $||\mbox{\boldmath$\bf L$}_{i}^{*\omega{\rm
MTM}[n](k+1)}-\mbox{\boldmath$\bf L$}_{i}^{*\omega{\rm MTM}[n](k)}||\leq{\rm
TOL}$ ($\exists k$) sufficiently close to zero and the iteration loop is
closed by again predicting a new state. If the
$\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}^{[n+1](k+1)}$
slowly converges, the increment
$\Delta\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}^{[n+1]}$
(4.45) should be decresed, e.g.,
$\Delta\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}^{[n+1]}:=\Delta\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}^{[n+1]}/2$.
This way an implicit iterative algorithm is set up to evaluate the stress
response to the average strain increment prescribed by Eq. (4.45).

We considered in Subsection 4.7 only one type of nonlinear peridynamics,
precisely nonlinear elasticity. However, a more detailed analyses of nonlinear
incremental methods for both different physical phenomena (e.g., diffusion,
viscosity, thermoelastoplasticity, damage accumulation, debonding, plastic
localization, wave propagation) and different random and periodic structures
(e.g., polycrystals, fiber networks, hybrid structures, foam materials) are
beyond the scope of the current study and deserve further investigations.

##  5 Estimation of field fluctuations and effective energy-based criteria

Although the effective properties of CM is the main goal of micromechanics, it
is also important to provide insight into the statistical distribution of
local stresses in each phase. It is not a problem if we have a direct
numerical simulation (DNS) as for periodic structure CMs (see, e.g., [120],
[204]). However, consideration of random structure CMs with nonlinear
constitutive behavior, including plasticity, damage, and creep, requires an
acceptance of some additional assumptions. The popular hypothesis of the mean-
field method takes into account only the average stresses in the phases
$\mbox{$\langle$}\mbox{\boldmath$\sigma$}\mbox{$\rangle$}_{i}$, which are used
in most secant methods (see, e.g., [29], [270]). However, the second
statistical moment of stress
$\mbox{$\langle$}\mbox{\boldmath$\sigma$}\otimes\mbox{\boldmath$\sigma$}\mbox{$\rangle$}_{i}(\mbox{\boldmath$\bf
x$})$ is each phase $\mbox{\boldmath$\bf x$}\in v_{i}\subset v^{(k)}\
(k=0,1,\ldots,n)$ is directly used in the analyses of a wide class of
nonlinear phenomena (strength, plasticity, and fracture) and some critical
effects (e.g., dielectric breakdown, see [80], [244]. So, in the strength
theory, the Tsa-Wu criterion has a tensor polynomial (see, e.g., [273])
failure surface depending on stresses. The onset of yielding in some cases of
elastic-plastic materials has also a quadratic yield condition (see, e.g.,
[27]) generalizing some classical criteria, as, for example, the von Mises
equivalent stress criterion. Stress potential in power-law creep is a
nonlinear function of equivalent stresses (see, e.g., [43], [60], [234]). In
linear fracture theory, the energy release rate for penny-shape crack has a
quadratic dependence on the stress intensity factors (see [242]) which, in
turn, are the linear functions of the stresses. One of the main achievements
in local nonlinear micromechanics (see [234]) is a suitable approximation in
the context of the variational principles permitting determining bounds for
nonlinear composites, directly from the corresponding bounds for linear
comparison composites with the same microstructures as the nonlinear
composites. It was also concluded (see [267]) that the variational method can
be interpreted as a secant method based on the estimations of the second
moment of stresses.

Conceptually, more general nonlinear models were proposed in peridynamics. So,
classical von Mises plasticity theory is generalized by Silling [260] by the
introduction of a scale-valued yield surface as a nonlinear function of force
state. In particular, the yield surface can be a quadratic function of the
norm expressed through a dot product of the force state (analog of an
equivalent Mises stress). In a peridynamic elastoplastic model proposed, a
deformation state is decomposed to isotropic and deviatoric parts whereas a
scale force state is decomposed on co-isotropic and co-deviatoric parts
allowing the condition of plastic incompressibility to be enforced. Mitchell
[212] developed this model by integrating the flow rule tracking along the
lines of the classical theories of rate-independent $J_{2}$ plasticity. The
author proposed a useful non-local yield criterion in terms of 3D force states
equivalent to $J_{2}$ plasticity depending on the yield stress and horizon for
the material. Madenci and Oterkus [198] presented constitutive relations for
plastic deformation according to von Mises yield criteria with linear
isotropic hardening. The expressions for the yield function and the rule of
incremental plastic stretch are derived in terms of the horizon, shear
modulus, force density, and hardening parameter of the medium. The yield
surface is constructed using the relationship between the effective stress and
equivalent plastic stretch. This model was later extended in [232] (see also
[9]) where the von Mises yield criteria are exploited for a description of
plastic yielding whereas the equivalent plastic stretch is utilized as an
internal variable for the general form of the mixture of isotropic and
kinematic hardening. One [220] formulated two rate-independent yield functions
equivalent to $J_{2}$ plasticity with associated flow rules using the force
states and the strain energy density for 2D cases.

Another area of nonlinear phenomena is the damage and failure where a bond is
being “broken” or “intact” following certain failure criteria. Damage modeling
initiated by the seminal work by Silling and Askari [259] (see also [197]) was
determined by using critical stretch criteria of bonds. The potential role of
employing a nonlinear peridynamic kernel in predicting the onset of fractures
has been explored in [96] (see also [84], [85]). It was demonstrated that the
presence of nonlinearity in the relative displacement can reproduce inherently
the high potential energy peaks that indicate the formation of cracks. Another
direction of damage modeling is associated with energy-based criteria such as,
e.g., the critical bond strain energy density failure criteria [124], [268],
[269], [279], [299]. The last model was applied in energy-based peridynamic
approach [224] for fatigue cracking proposed. Here, the energy release rate
for interaction between material points is expressed through the
micropotentials of the interaction between these points. The definition of the
critical energy release rate range and the energy-based peridynamic failure
model for both phases crack initiation and crack growth phases are introduced.
Pointwise energy-based bond breaking model states (see [150]) that a bond is
broken and no longer carries load if the strain energy density at either of
its endpoints exceeds some critical local strain energy density $W_{c}$ when
the point loses all of its bonds and is completely detached from the body.

In the nonlinear problems of the LM considered, the central deal is an
estimation of the second moment of phase fields
($\mbox{\boldmath$\pi$}=\mbox{\boldmath$\varepsilon$},\
\mbox{\boldmath$\sigma$}$; $\mbox{\boldmath$\bf x$}\in v^{(i)}$)

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\langle\mbox{\boldmath$\pi$}\otimes\mbox{\boldmath$\pi$}\rangle_{i}({\bf
x})=\langle\mbox{\boldmath$\pi$}\rangle_{i}({\bf
x})\otimes\langle\mbox{\boldmath$\pi$}\rangle_{i}({\bf x})+\Delta^{\pi
2}_{i}(\mbox{\boldmath$\bf x$}).$ (5.1)

where $\Delta^{\pi 2}_{i}(\mbox{\boldmath$\bf x$})$ is called the dispersion
(or the fluctuation field). There are two methods for estimations of
$\langle\mbox{\boldmath$\pi$}\otimes\mbox{\boldmath$\pi$}\rangle_{i}({\bf
x})$: the exact differential method (see pp. 419, 423 in [43], and pp. 252,
253 in [60]), and the method of integral equations (see pp. 427-431 in [43],
and pp. 446, 456-475, 532 in [60]). Both methods lead to close numerical
results. It was also demonstrated that engineering approach
$\langle\mbox{\boldmath$\pi$}\otimes\mbox{\boldmath$\pi$}\rangle_{i}({\bf
x})=\langle\mbox{\boldmath$\pi$}\rangle_{i}({\bf
x})\otimes\langle\mbox{\boldmath$\pi$}\rangle_{i}({\bf x})$ can lead to an
infinite error at the hydrostatic loading.

However, it is precisely the method of integral equations in LM that enables
one to generalize it to the PM. Namely, the following new approximations for
the second moment can be obtained by taking only the binary interaction of
inclusions into account:

$\displaystyle\langle\overline{}\mbox{\boldmath$\vartheta$}\otimes\overline{}\mbox{\boldmath$\vartheta$}\rangle_{i}({\bf
x})$ $\displaystyle=$
$\displaystyle\langle\overline{}\mbox{\boldmath$\vartheta$}\rangle_{i}({\bf
x})\otimes\langle\overline{}\mbox{\boldmath$\vartheta$}\rangle_{i}({\bf
x})+\int\\!\\!\mbox{$\langle$}\mbox{\boldmath$\cal
L$}_{p}^{\theta\zeta}(\mbox{\boldmath$\bf
x$},\overline{\mbox{\boldmath$\zeta$}})|;v_{i},\mbox{\boldmath$\bf
x$}_{i}\mbox{$\rangle$}$ (5.2) $\displaystyle\otimes$
$\displaystyle\mbox{$\langle$}\\!\mbox{\boldmath$\cal
L$}_{p}^{\theta\zeta}(\mbox{\boldmath$\bf
x$},\overline{\mbox{\boldmath$\zeta$}})|;v_{i},{\bf
x}_{i}\rangle\varphi(v_{p},{\bf x}_{p}\mid;v_{i},{\bf x}_{i})d{\bf x}_{p}.$

Using the known properties of the second statistical moment, we obtain the
relations for the average second moments of the fields,
$\langle\mbox{\boldmath$\vartheta$}\otimes\mbox{\boldmath$\vartheta$}\rangle_{i}({\bf
x})$ in the component $v_{i}$ corresponding to Eqs. (5.2)

$\displaystyle\\!\\!\langle\mbox{\boldmath$\vartheta$}\otimes\mbox{\boldmath$\vartheta$}\rangle_{i}({\bf
x})$ $\displaystyle=$
$\displaystyle\langle\mbox{\boldmath$\vartheta$}\rangle_{i}({\bf
x})\otimes\langle\mbox{\boldmath$\vartheta$}\rangle_{i}({\bf
x})+\\!\\!\int\\!\\!\mbox{\boldmath$\cal
A$}^{\theta\theta}_{i}*\mbox{$\langle$}\mbox{\boldmath$\cal
L$}_{p}^{\theta\zeta}(\mbox{\boldmath$\bf
x$},\overline{\mbox{\boldmath$\zeta$}})|;v_{i},\mbox{\boldmath$\bf
x$}_{i}\mbox{$\rangle$}$ (5.3) $\displaystyle\otimes$
$\displaystyle\mbox{\boldmath$\cal
A$}^{\theta\theta}_{i}*\mbox{$\langle$}\\!\mbox{\boldmath$\cal
L$}_{p}^{\theta\zeta}(\mbox{\boldmath$\bf
x$},\overline{\mbox{\boldmath$\zeta$}})|;v_{i},{\bf
x}_{i}\rangle\varphi(v_{p},{\bf x}_{p}\mid;v_{i},{\bf x}_{i})d{\bf x}_{p}.$

Comment 5.1. New equations (5.2) and (5.3) are obtained for nonlinear problems
of PM with nonlinear operators $\mbox{\boldmath$\cal A$}^{\theta\theta}_{i}$
(3.232) and $\mbox{\boldmath$\cal L$}_{p}^{\theta\zeta}(\mbox{\boldmath$\bf
x$},\overline{\mbox{\boldmath$\zeta$}})$ (3.231). These equations are reduced
to the corresponding equations for the linear bond-based model [65] and linear
state-based model [67] as well as to the equations obtained before for the
linear LM (in the framework of both the first background [43] and second
background [60] of micromechanics). Thus, we have the third (see Comments
3.3.2 and 4.4.) rare case when the methods developed in PM for the field
second moments (nonlinear problem) initiate the new methods of LM (a
generalization of linear analyses to nonlinear ones) unknown before. In so
doing, linearity of matrix properties as in Comment 4.4 is not required.
Moreover, Eqs. (5.2) and(5.3) including double-point $\varphi(v_{p},{\bf
x}_{p}|;v_{i},{\bf x}_{i})$ conditional probability densities can be
generalized (analogously to [67]) to the corresponding exact representations
containing triple-point $\varphi(v_{p},{\bf x}_{p},v_{q},{\bf
x}_{q}|;v_{i},{\bf x}_{i})$ densities.

Simplifications of Eq. (5.3) with subsequent applications to obtaining of
effective energy-based criteria can be performed for linear operators
$\mbox{\boldmath$\cal A$}^{\theta\theta}_{i}$ (3.232) and
$\mbox{\boldmath$\cal L$}_{p}^{\theta\zeta}(\mbox{\boldmath$\bf
x$},\overline{\mbox{\boldmath$\zeta$}})$ (3.231) in the framework of the EFH
$\bf H$1a (4.1):

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\langle\mbox{\boldmath$\eta$}\otimes\mbox{\boldmath$\eta$}\rangle_{i}({\mbox{\boldmath$\bf
z$}})$ $\displaystyle=$
$\displaystyle\langle\mbox{\boldmath$\eta$}\rangle_{i}({\bf
z})\otimes\langle\mbox{\boldmath$\eta$}\rangle_{i}({\bf z})+\Delta^{\eta
2}_{i}(\mbox{\boldmath$\bf z$}),$ (5.4)

where

$\displaystyle\\!\\!\\!\\!\\!\\!\Delta^{\eta 2}_{i}(\mbox{\boldmath$\bf z$})$
$\displaystyle=$
$\displaystyle\int[\mbox{$\mathbb{A}$}_{i,j}(\mbox{\boldmath$\bf
z$})\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}]\otimes[\mbox{$\mathbb{A}$}_{i,j}(\mbox{\boldmath$\bf
z$})\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}]\varphi(v_{j},{\bf
x}_{j}\mid;v_{i},{\bf x}_{i})~{}d\mbox{\boldmath$\bf x$}_{j},$
$\displaystyle\\!\\!\\!\\!\\!\\!\mbox{\boldmath$\bf
A$}_{i,j}(\mbox{\boldmath$\bf z$})$ $\displaystyle=$
$\displaystyle\mbox{\boldmath$\bf
A$}^{\eta\varepsilon}_{i}(\mbox{\boldmath$\bf z$})[\mbox{\boldmath$\bf
J$}^{I\varepsilon\varepsilon}_{i,j}\mbox{\boldmath$\bf
D$}^{0}_{i}+\mbox{\boldmath$\bf
J$}^{J\varepsilon\varepsilon}_{i,j}(\mbox{\boldmath$\bf
y$})\mbox{\boldmath$\bf D$}^{0}_{j}].$ (5.5)

Here $\mbox{\boldmath$\bf A$}^{\eta\varepsilon}_{i}(\mbox{\boldmath$\bf
z$})=\underline{\mbox{\boldmath$\bf
A$}}_{i}^{u\varepsilon}[\mbox{\boldmath$\bf
x$}]\mbox{$\langle$}\mbox{\boldmath$\xi$}\mbox{$\rangle$}$
($\mbox{\boldmath$\xi$}=\hat{}\mbox{\boldmath$\bf x$}-\mbox{\boldmath$\bf
x$})$ whereas $\mbox{\boldmath$\bf D$}^{0}_{i}$ is the effective field
concentration factor
$\mbox{$\langle$}\overline{\mbox{\boldmath$\varepsilon$}}\mbox{$\rangle$}_{i}=\mbox{\boldmath$\bf
D$}^{0}_{i}\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}$ (see
for details [67]).

Equation (5.4) is a “elementary block” for the construction of many effective
energy-based criteria. So, in the energy-based peridynamic approach for
fatigue cracking, one [224] estimated the energy release rate for interaction
between material points through the micropotentials of the interaction between
these points. The state-based version of the effective energy-based criterion
is the following $(\mbox{\boldmath$\xi$},\mbox{\boldmath$\zeta$}\in{\cal
H}_{\bf x}\subset v_{i})$

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\langle\underline{\mbox{\boldmath$\bf
U$}}^{\top}\underline{\mbox{$\mathbb{K}$}}\underline{\mbox{\boldmath$\bf
U$}}\rangle_{i}(\mbox{\boldmath$\xi$},\mbox{\boldmath$\zeta$})$
$\displaystyle=$ $\displaystyle\langle\underline{\mbox{\boldmath$\bf
U$}}^{\top}\mbox{$\rangle$}_{i}\underline{\mbox{$\mathbb{K}$}}\mbox{$\langle$}\underline{\mbox{\boldmath$\bf
U$}}\rangle_{i}(\mbox{\boldmath$\xi$},\mbox{\boldmath$\zeta$})+\Delta^{\rm
UKU}(\mbox{\boldmath$\xi$},\mbox{\boldmath$\zeta$}),$ (5.6)
$\displaystyle\Delta^{\rm UKU}(\mbox{\boldmath$\xi$},\mbox{\boldmath$\zeta$})$
$\displaystyle=$
$\displaystyle\int[\mbox{$\mathbb{A}$}_{ij}(\mbox{\boldmath$\xi$})\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}]^{\top}\underline{\mbox{$\mathbb{K}$}}(\mbox{\boldmath$\xi$},\mbox{\boldmath$\zeta$})[\mbox{$\mathbb{A}$}_{ij}(\mbox{\boldmath$\zeta$})\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}]\varphi(v_{j},{\bf
x}_{j}\mid;v_{i},{\bf x}_{i})~{}d\mbox{\boldmath$\bf x$}_{j},$ (5.7)

where the statistical average of the displacement state is defined as
$\mbox{$\langle$}\underline{\mbox{\boldmath$\bf
U$}}\rangle_{i}[\mbox{\boldmath$\bf
x$}]\mbox{$\langle$}\mbox{\boldmath$\xi$}\mbox{$\rangle$}=\underline{\mbox{\boldmath$\bf
A$}}^{u\epsilon*}_{i}[\mbox{\boldmath$\bf
x$}]\mbox{$\langle$}\mbox{\boldmath$\xi$}\mbox{$\rangle$}\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}=\mbox{\boldmath$\bf
A$}^{\eta\epsilon*}_{i}(\mbox{\boldmath$\bf
z$})\mbox{$\langle$}\mbox{\boldmath$\xi$}\mbox{$\rangle$}\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}$
(4.17). The energy-based failure model [150] is based on an evaluation of the
local strain energy density for bond-based peridynamics at each point. Then
the effective criterion for the state-based model can be obtained

$\displaystyle\\!\\!\\!4\mbox{$\langle$}{\cal
W}\mbox{$\rangle$}_{i}(\mbox{\boldmath$\bf x$})$ $\displaystyle\\!\\!=\\!\\!$
$\displaystyle\langle\underline{\mbox{\boldmath$\bf
U$}}^{\top}\bullet\underline{\mbox{$\mathbb{K}$}}\bullet\underline{\mbox{\boldmath$\bf
U$}}\rangle_{i}({\bf x})=\langle\underline{\mbox{\boldmath$\bf
U$}}^{\top}\mbox{$\rangle$}_{i}\bullet\underline{\mbox{$\mathbb{K}$}}\bullet\mbox{$\langle$}\underline{\mbox{\boldmath$\bf
U$}}\rangle_{i}({\bf x})+\Delta^{\rm UKU}(\mbox{\boldmath$\bf x$}),$ (5.8)
$\displaystyle\Delta^{\rm UKU}(\mbox{\boldmath$\bf x$})$ $\displaystyle=$
$\displaystyle\\!\int\\!\\!\Big{\langle}[\mbox{$\mathbb{A}$}_{ij}(\mbox{\boldmath$\xi$})\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}]^{\top}\big{\langle}\underline{\mbox{$\mathbb{K}$}}(\mbox{\boldmath$\xi$},\mbox{\boldmath$\zeta$})[\mbox{$\mathbb{A}$}_{ij}(\mbox{\boldmath$\zeta$})\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}]\big{\rangle}^{\cal
H_{\bf x}}\Big{\rangle}^{\cal H_{\bf x}}\\!\\!(\mbox{\boldmath$\bf x$})$ (5.9)
$\displaystyle\times$ $\displaystyle\varphi(v_{j},{\bf x}_{j}\mid;v_{i},{\bf
x}_{i})~{}d\mbox{\boldmath$\bf x$}_{j}$

averaging over the horizon region, ${\cal H}_{\bf x}\subset v_{i}$
($\mbox{\boldmath$\bf y$}\in v_{j}$).

Seemingly, a second moment of the force vector state
$\underline{\mbox{\boldmath$\bf
T$}}=\underline{\mbox{$\mathbb{K}$}}\bullet\underline{\mbox{\boldmath$\bf
U$}}$ (2.19) can be obtained
$(\mbox{\boldmath$\xi$},\mbox{\boldmath$\zeta$}\in{\cal H}_{\bf x}\subset
v_{i}$)

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\mbox{$\langle$}\underline{\mbox{\boldmath$\bf
T$}}\otimes\underline{\mbox{\boldmath$\bf T$}}\mbox{$\rangle$}_{i}^{\cal
H_{\bf x}}(\mbox{\boldmath$\bf x$})$ $\displaystyle=$
$\displaystyle\Big{\langle}(\underline{\mbox{$\mathbb{K}$}}\bullet\mbox{$\langle$}\underline{\mbox{\boldmath$\bf
U$}}\mbox{$\rangle$}_{i})\otimes(\underline{\mbox{$\mathbb{K}$}}\bullet\mbox{$\langle$}\underline{\mbox{\boldmath$\bf
U$}}\mbox{$\rangle$}_{i})\Big{\rangle}^{\cal H_{\bf x}}(\mbox{\boldmath$\bf
x$})+\Delta^{\rm T2}(\mbox{\boldmath$\bf x$}),$ (5.10)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\Delta^{\rm
T2}(\mbox{\boldmath$\bf x$})$ $\displaystyle=$
$\displaystyle\int\\!\\!\Big{\langle}\big{\langle}\underline{\mbox{$\mathbb{K}$}}(\mbox{\boldmath$\xi$},\mbox{\boldmath$\zeta$})[\mbox{$\mathbb{A}$}_{ij}(\mbox{\boldmath$\zeta$})\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}]\big{\rangle}^{\cal
H_{\bf
x}}\otimes\big{\langle}\underline{\mbox{$\mathbb{K}$}}(\mbox{\boldmath$\xi$},\mbox{\boldmath$\zeta$})$
(5.11) $\displaystyle\times$
$\displaystyle[\mbox{$\mathbb{A}$}_{ij}(\mbox{\boldmath$\zeta$})\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}]\big{\rangle}^{\cal
H_{\bf x}}\Big{\rangle}^{\cal H_{\bf x}}(\mbox{\boldmath$\bf
x$})\varphi(v_{j},{\bf x}_{j}\mid;v_{i},{\bf x}_{i})~{}d\mbox{\boldmath$\bf
x$}_{j}.$

After that, analogously to Eq. (5.8), we can find from Eq. (5.10)
$\mbox{$\langle$}\underline{\mbox{\boldmath$\bf
T$}}\cdot\underline{\mbox{\boldmath$\bf T$}}\mbox{$\rangle$}_{i}^{\cal H_{\bf
x}}(\mbox{\boldmath$\bf x$})$ ($\mbox{\boldmath$\bf x$}\in{\cal H}_{\bf
x}\subset v_{i}$). We can also estimate a particular case of the yield
function based on the criterium of the co-deviatoric force state:
$\psi(\mbox{\boldmath$\bf x$})=\big{\langle}\underline{t}^{\rm
d}[x]\mbox{$\langle$}\xi\mbox{$\rangle$}$ $\times\underline{t}^{\rm
d}[x]\mbox{$\langle$}\xi\mbox{$\rangle$}\rangle^{\cal H_{\bf
x}}(\mbox{\boldmath$\bf x$})/2$ (which is equivalent to the second invariant
of the deviatoric stress tensor in classical local theory), where
$\underline{t}^{\rm d}[x]\mbox{$\langle$}\xi\mbox{$\rangle$}$ is the co-
deviatoric part of the scalar force state $\underline{t}$ (see for details
[212], [260]).

A failure criterion for breaking the physical bonds between material points
based on the energy dissipation ${\cal W}_{\bf\xi}(\mbox{\boldmath$\bf x$})$
(see [124], [268], [269]) is a function of the stress state
$\underline{\mbox{\boldmath$\bf
T$}}=\underline{\mbox{$\mathbb{K}$}}\bullet\underline{\mbox{\boldmath$\bf
U$}}$ (2.19) and the displacement state $\underline{\mbox{\boldmath$\bf U$}}$
between two points ($\mbox{\boldmath$\xi$}\in{\cal H}_{\bf x}\subset v_{i})$

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!{\cal
W}_{\bf\xi}(\mbox{\boldmath$\bf x$})\\!$ $\displaystyle=$
$\displaystyle\\!\big{\langle}(\underline{\mbox{$\mathbb{K}$}}\bullet\underline{\mbox{\boldmath$\bf
U$}})\rangle_{i}\langle\underline{\mbox{\boldmath$\bf
U$}}\big{\rangle}_{i}(\mbox{\boldmath$\bf
x$})=(\underline{\mbox{$\mathbb{K}$}}^{(i)}\bullet\langle\underline{\mbox{\boldmath$\bf
U$}}\rangle_{i})\langle\underline{\mbox{\boldmath$\bf
U$}}\rangle_{i}(\mbox{\boldmath$\bf x$})+\Delta^{\rm
TU2}(\mbox{\boldmath$\xi$},\mbox{\boldmath$\bf x$}),$ (5.12)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\Delta^{\rm
TU2}(\mbox{\boldmath$\xi$},\mbox{\boldmath$\bf x$})\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\int\\!\big{\langle}\underline{\mbox{$\mathbb{K}$}}(\mbox{\boldmath$\xi$},\mbox{\boldmath$\zeta$})[\mbox{$\mathbb{A}$}_{ij}(\mbox{\boldmath$\zeta$})\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}]\big{\rangle}^{\cal
H_{\bf
x}}[\mbox{$\mathbb{A}$}_{ij}(\mbox{\boldmath$\xi$})\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}](\mbox{\boldmath$\bf
x$})\varphi(v_{j},\\!{\bf x}_{j}\\!\mid;v_{i},{\bf x}_{i})d\mbox{\boldmath$\bf
x$}_{j}.$ (5.13)

Equations (5.6)-(5.13) for the linear operators $\mbox{\boldmath$\cal
A$}^{\theta\theta}_{i}$ (3.232) and $\mbox{\boldmath$\cal
L$}_{p}^{\theta\zeta}(\mbox{\boldmath$\bf
x$},\overline{\mbox{\boldmath$\zeta$}})$ (3.231). are reduced to the
corresponding equations for both the linear state-based [67] and linear bond-
based model [65] models

Comment 5.2. It is interesting that all effective energy-based criteria
considered (5.6)-(5.13) are constructed by the use of the same “elementary
block” of the second moment of the relative displacement
$\mbox{$\langle$}\mbox{\boldmath$\eta$}\otimes\mbox{\boldmath$\eta$}\mbox{$\rangle$}_{i}(\mbox{\boldmath$\bf
z$})$ (5.4) (or, that is the same,
$\mbox{$\langle$}\underline{\mbox{\boldmath$\bf U$}}[\mbox{\boldmath$\bf
x$}]\mbox{$\langle$}\mbox{\boldmath$\xi$}\mbox{$\rangle$}\otimes\underline{\mbox{\boldmath$\bf
U$}}[\mbox{\boldmath$\bf
x$}]\mbox{$\langle$}\mbox{\boldmath$\xi$}\mbox{$\rangle$}\mbox{$\rangle$}_{i}(\mbox{\boldmath$\bf
x$})$ ($\mbox{\boldmath$\bf x$}\in v_{i}$, $\mbox{\boldmath$\xi$}\in{\cal
H}_{\bf x}\subset v_{i}$)). Moreover, these effective energy-based criteria
are defined on either one bond $\mbox{\boldmath$\xi$}\in{\cal H}_{\bf
x}\subset v_{i}$, or two bonds
$\mbox{\boldmath$\xi$},\mbox{\boldmath$\zeta$}\in{\cal H}_{\bf x}\subset
v_{i}$ with the field fluctuations, either $\Delta^{\rm
UKU}(\mbox{\boldmath$\xi$},\mbox{\boldmath$\zeta$}),\Delta^{\rm
UKU}(\mbox{\boldmath$\bf x$})$, $\Delta^{\rm T2}(\mbox{\boldmath$\bf x$})$, or
$\Delta^{\rm TU2}(\mbox{\boldmath$\xi$},\mbox{\boldmath$\bf x$})$.

Comment 5.3. The scheme of estimations of different effective energy-based
criteria by the use of one “elementary block”
$\langle\mbox{\boldmath$\eta$}\otimes\mbox{\boldmath$\eta$}\rangle_{i}({\mbox{\boldmath$\bf
z$}})$ (5.4) is presented in Fig. 6.

Fig. 6: Estimation of the effective energy-based criteria by the use of the
elementary block
$\mbox{$\langle$}\mbox{\boldmath$\eta$}\otimes\mbox{\boldmath$\eta$}\mbox{$\rangle$}_{i}(\mbox{\boldmath$\bf
z$})$

##  6 Fast Fourier transform in peridynamic micromechanics

###  6.1 Fast Fourier transform

We consider a unit cell $\Omega_{00}={\rm\Pi}_{\alpha}[0,l^{\Omega}_{\alpha}]$
($\alpha=1,\ldots,d)$ discretized with the uniform grid of
$N_{1}\times\ldots\times N_{d}$ nodes

$\displaystyle\mbox{\boldmath$\bf x$}^{\bf k}_{\bf N}$ $\displaystyle=$
$\displaystyle\mbox{\boldmath$\bf k$}\cdot\Delta\mbox{\boldmath$\bf
h$}\equiv[k_{1}\Delta h_{1},\ldots,k_{d}\Delta h_{d}]^{\top},$ $\displaystyle\
\ \mbox{\boldmath$\bf k$}\in Z^{d}_{\bf N}$ $\displaystyle:=$
$\displaystyle\Big{\\{}\mbox{\boldmath$\bf m$}\in Z^{d}\Big{|}0\leq
m_{\alpha}\leq N_{\alpha}-1\Big{\\}},$ (6.1)

hereafter $Z^{d}$ denotes the set of $d-$tuples of integers, $N_{\alpha}$
($\alpha=1,\ldots,d)$ are odd, $\mbox{\boldmath$\bf
N$}:=(N_{1},\ldots,N_{d})^{\top}$, $|\mbox{\boldmath$\bf
N$}|={\rm\Pi}_{\alpha}N_{\alpha}$, and $\Delta
h_{\alpha}=l^{\Omega}_{\alpha}/N_{\alpha}$ are the steps of grid in the
direction $\alpha$ with the unite vector $\mbox{\boldmath$\bf e$}_{\alpha}$.
In this section, we use the notation
$\Omega_{00}={\rm\Pi}_{\alpha}[0,l^{\Omega}_{\alpha}]$ popular in the DFT (and
the corresponding indexation) instead of the notation
$\Omega_{00}=\Pi_{\alpha}[-l^{\Omega}_{\alpha},l^{\Omega}_{\alpha}]$ exploited
in PM.

The discrete Fourier transform $\widehat{\mbox{\boldmath$\bf f$}}=\mbox{$\cal
F$}_{d}(\mbox{\boldmath$\bf f$})$ (DFT) and the inverse DFT (iDFT)
$\mbox{\boldmath$\bf f$}=\mbox{$\cal F$}_{d}^{-1}(\widehat{\mbox{\boldmath$\bf
f$}})$ for the discrete periodic functions are defined as (see [37], [38],
[203])

$\displaystyle\widehat{\mbox{\boldmath$\bf f$}}_{\bf N}(\mbox{\boldmath$\bf
k$})=\sum_{{\bf n}\in Z^{d}_{\bf N}}\mbox{\boldmath$\bf f$}_{\bf
N}(\mbox{\boldmath$\bf n$})\omega^{-{\bf k}{\bf n}}_{\bf N},\ \ \
\mbox{\boldmath$\bf f$}_{\bf N}(\mbox{\boldmath$\bf
n$})=\frac{1}{|\mbox{\boldmath$\bf N$}|}\sum_{{\bf k}\in Z^{d}_{\bf
N}}\widehat{\mbox{\boldmath$\bf f$}}_{\bf N}(\bf k)\omega^{{\bf k}{\bf
n}}_{\bf N}$ (6.2)

for $\bf N$-point vectors (called also sequences) $\mbox{\boldmath$\bf
k$},\mbox{\boldmath$\bf n$}\in Z^{d}_{\bf N}$; hereafter $\mbox{\boldmath$\bf
f$}(\mbox{\boldmath$\bf n$})$ and $\widehat{\mbox{\boldmath$\bf
f$}}(\mbox{\boldmath$\bf k$})$ are the sequences of $\bf N$ possibly complex-
valued numbers, and $\omega^{{\bf k}{\bf n}}_{\bf N}={\rm exp}[2\pi
i\sum_{\alpha}(k_{\alpha}n_{\alpha})/N_{\alpha}]$, where $\alpha=1,\ldots,d$
and $i=\sqrt{-1}$.

The properties of the DFT that are most applied in subsequent presentations
are

$\displaystyle\mbox{$\cal F$}_{d}(\nabla\mbox{\boldmath$\bf
f$})=i\mbox{\boldmath$\bf k$}\widehat{}\mbox{\boldmath$\bf f$},\ \ \
\mbox{$\cal F$}_{d}(\mbox{\boldmath$\bf f$}\odot\mbox{\boldmath$\bf g$})_{\bf
N}(\mbox{\boldmath$\bf k$})=\widehat{\mbox{\boldmath$\bf f$}}_{\bf
N}(\mbox{\boldmath$\bf k$})\cdot\widehat{\mbox{\boldmath$\bf g$}}_{\bf
N}(\mbox{\boldmath$\bf k$}),$ (6.3)

where

$\displaystyle(\mbox{\boldmath$\bf f$}\odot\mbox{\boldmath$\bf
g$})(\mbox{\boldmath$\bf x$}^{\bf n}_{\bf N}):=\sum_{\bf j}\mbox{\boldmath$\bf
f$}(\mbox{\boldmath$\bf x$}^{\bf n-\bf j}_{\bf N})\cdot\mbox{\boldmath$\bf
g$}(\mbox{\boldmath$\bf x$}^{\bf j}_{\bf N})$ (6.4)

is a circular (also called periodic) convolution. To obtain a nonperiodic (or
noncircular or aperiodic) convolution of two (one length $\bf L$ and another
length $\bf M$) sequences using DFT/iDFT approach, the two sequences
$\mbox{\boldmath$\bf f$}(\mbox{\boldmath$\bf x$}_{\bf L}^{\bf n})$ and
$\mbox{\boldmath$\bf g$}(\mbox{\boldmath$\bf x$}_{\bf M}^{\bf n})$ have to be
extended by adding zeros at the ends such as their length are
$\mbox{\boldmath$\bf N$}_{\rm e}\geq\mbox{\boldmath$\bf
L$}+\mbox{\boldmath$\bf M$}-{\bf 1}$ (This process is called zero padding).
So, for $L$-point sequence in $R^{1}$
$\\{g(n)\\}=\\{g_{0},g_{1},\ldots,g_{L-1}\\}$, the extended (zero-padded)
sequence is $\\{g_{\rm e}(n)\\}=\\{g_{0},g_{1},\ldots,g_{L-1},0,\ldots,0\\}$,

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\{g(n)\\}=\\{g_{0},g_{1},\ldots,g_{L-1}\\}\
\\!\longrightarrow\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!^{\rm ZP}\\!\ \ \ \\{g_{\rm
e}(n)\\}=\\{g_{0},g_{1},\ldots,g_{L-1},0,\ldots,0\\},$ (6.5)

where $\\{g_{\rm e}(n)\\}=0$ at $L\leq n\leq N_{\rm e}-1$. Without loss of
generality, we will assume in the current study that the length
$\mbox{\boldmath$\bf N$}_{\rm e}=\mbox{\boldmath$\bf L$}+\mbox{\boldmath$\bf
M$}$ is a power of 2 and $\mbox{\boldmath$\bf L$}=\mbox{\boldmath$\bf
M$}=\mbox{\boldmath$\bf N$}$. In such a case, e.g., Eqs. (6.21) and (6.5) can
be recast in the form

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\widehat{\mbox{\boldmath$\bf
f$}}_{\rm e{\bf N}}(\mbox{\boldmath$\bf k$})=\sum_{{\bf n}\in Z^{d}_{{\bf
N}_{\rm e}}}\mbox{\boldmath$\bf f$}_{\rm e{\bf N}}(\mbox{\boldmath$\bf
n$})\omega^{-{\bf k}{\bf n}}_{\bf N},\ \ (\mbox{\boldmath$\bf f$}_{\rm
e}\odot\mbox{\boldmath$\bf g$}_{\rm e})(\mbox{\boldmath$\bf x$}^{\bf n}_{{\bf
N}_{\rm e}}):=\sum_{\bf j}\mbox{\boldmath$\bf f$}_{\rm e}(\mbox{\boldmath$\bf
x$}^{\bf n-\bf j}_{{\bf N}_{\rm e}})\cdot\mbox{\boldmath$\bf g$}_{\rm e}^{\bf
j}(\mbox{\boldmath$\bf x$}^{\bf j}_{{\bf N}_{\rm e}})$ (6.6)

exploiting both the extended sequences $\mbox{\boldmath$\bf f$}_{\rm
e}(\mbox{\boldmath$\bf n$})$ and $\mbox{\boldmath$\bf g$}_{\rm
e}(\mbox{\boldmath$\bf n$})$, and the modified length $\mbox{\boldmath$\bf
N$}_{\rm e}$ depending on the lengths $\bf L$ and $\bf M$ of the initial
sequences $\mbox{\boldmath$\bf f$}(\mbox{\boldmath$\bf n$})$ and
$\mbox{\boldmath$\bf g$}(\mbox{\boldmath$\bf n$})$. In general, the DFT is a
very expensive algorithm $(O(N^{2}))$ requiring $N$ operations for each of the
$N$ Fourier components. However, this cost can be significantly reduced to the
computational cost of $O(N{\rm log}_{2}N)$ (where $N$ is usually taken as a
power of 2) as in the so-called Fast Fourier Transformation (FFT) algorithm
[87] (see for references [191]).

###  6.2 Fast Fourier transform method in peridynamic micromechanics

Expansion of the FFT methods into the LM was initiated by Moulinec and Suquet
[218], [219]. They performed the DFT (6.8) of the implicit Lippmann-Schwinger
(L-S) equation (6.7) (compare with Eq. (3.3))

$\displaystyle\mbox{\boldmath$\varepsilon$}(\mbox{\boldmath$\bf x$})$
$\displaystyle=$
$\displaystyle\mbox{\boldmath$\varepsilon$}^{w_{\Gamma}}+\mbox{\boldmath$\bf
U$}^{(0)}*^{L}\\!\mbox{\boldmath$\tau$}(\mbox{\boldmath$\bf x$}),$ (6.7)
$\displaystyle\widehat{}\mbox{\boldmath$\varepsilon$}(\mbox{\boldmath$\bf
k$})$ $\displaystyle=$ $\displaystyle\widehat{}\mbox{\boldmath$\bf
U$}^{(0)}(\mbox{\boldmath$\bf
k$})^{L}\\!\widehat{}\mbox{\boldmath$\tau$}(\mbox{\boldmath$\bf k$})\ \ \
(\mbox{\boldmath$\bf k$}\not={\bf
0}),\widehat{\mbox{\boldmath$\varepsilon$}}^{w_{\Gamma}}({\bf
0})={\mbox{\boldmath$\varepsilon$}}^{w_{\Gamma}},$ (6.8)

where $*$ stands for the convolution operation. The origin FFT based basic
scheme proposed [218], [219] was solved by the fixed-point iterative method
(Picard’s iterations).

It is interesting that Eq. (6.7) (with accuracy to notations) looks like
linear state-based peridynamic Eq. (2.21) and (2.22). Hence, there is the
temptation to use the same method to obtain the DFT of Eq. (2.22) (like, Eq.
(6.8)). Unfortunately, we have a few obstacles to the realization of this
temptation: O1) a displacement field $\mbox{\boldmath$\bf
u$}(\mbox{\boldmath$\bf x$})$ (2.22) is not $\Omega_{00}$-periodic; O2) the
kernel $\mbox{\boldmath$\bf C$}(\mbox{\boldmath$\bf x$},\mbox{\boldmath$\bf
q$})$ (2.22) is not a convolution for inhomogeneous media
($\mbox{\boldmath$\bf C$}(\mbox{\boldmath$\bf x$},\mbox{\boldmath$\bf
q$})\not=\mbox{\boldmath$\bf C$}(\mbox{\boldmath$\bf x$}-\mbox{\boldmath$\bf
q$})$), and O3) a convolution (even if we have it) is not a circular
convolution (6.4). We demonstrate how to overcome these obstacles.

Namely, in the so-called decomposition approach based on ‘extraction’ of a
field corresponding to the homogeneous displacement loading, we split the
displacement field

$\displaystyle\mbox{\boldmath$\bf u$}(\mbox{\boldmath$\bf
x$})=\mbox{\boldmath$\bf u$}_{0}(\mbox{\boldmath$\bf x$})+\mbox{\boldmath$\bf
u$}_{1}(\mbox{\boldmath$\bf x$}),\ \ \mbox{\boldmath$\bf
u$}_{0}(\mbox{\boldmath$\bf
x$}):=\mbox{\boldmath$\varepsilon$}^{w_{\Gamma}}\mbox{\boldmath$\bf x$},$
(6.9)

into its fluctuating $\Omega_{00}$-periodic part $\mbox{\boldmath$\bf
u$}_{1}(\mbox{\boldmath$\bf x$})$ with zero average and a linearly growing
part $\mbox{\boldmath$\bf u$}_{0}(\mbox{\boldmath$\bf x$})$ corresponding for
the macroscopic strain
$\mbox{\boldmath$\varepsilon$}^{w\Gamma}(\mbox{\boldmath$\bf x$})$. In
additionally to the convolution kernel $\mbox{\boldmath$\bf C$}^{(0)}$ we
introduce two convolution kernels

$\displaystyle\mbox{\boldmath$\bf C$}^{0}_{1}(\mbox{\boldmath$\bf
x$}-\mbox{\boldmath$\bf q$})$ $\displaystyle=$
$\displaystyle\mbox{\boldmath$\bf C$}^{(1)}(\mbox{\boldmath$\bf
x$}-\\!\mbox{\boldmath$\bf q$})-\mbox{\boldmath$\bf
C$}^{(0)}(\mbox{\boldmath$\bf x$}-\\!\mbox{\boldmath$\bf q$}),$
$\displaystyle\mbox{\boldmath$\bf C$}^{0}_{1\epsilon}(\mbox{\boldmath$\bf
x$}-\\!\mbox{\boldmath$\bf q$})$ $\displaystyle=$
$\displaystyle\mbox{\boldmath$\bf C$}^{0}_{1}(\mbox{\boldmath$\bf
x$}-\mbox{\boldmath$\bf
q$})\mbox{\boldmath$\varepsilon$}^{w_{\Gamma}}(\mbox{\boldmath$\bf
q$}-\mbox{\boldmath$\bf x$})$ (6.10)

with the same horison $l_{\delta}$ as the horison of $\mbox{\boldmath$\bf
C$}^{(0)}(\mbox{\boldmath$\bf x$}-\mbox{\boldmath$\bf q$})$. It allows us to
decompose the equilibrium equation as

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!{\mbox{\boldmath$\cal
L$}}^{(0)}(\mbox{\boldmath$\bf C$}^{(0)},\mbox{\boldmath$\bf
u$}_{1})(\mbox{\boldmath$\bf x$})\\!=\\!{\mbox{\boldmath$\cal
L$}}_{1}(\mbox{\boldmath$\bf C$}_{1},\mbox{\boldmath$\bf
u$}_{1})(\mbox{\boldmath$\bf x$})-\mbox{\boldmath$\bf
b$}^{w\Gamma}(\mbox{\boldmath$\bf x$}),\ \mbox{\boldmath$\bf
b$}^{w\Gamma}(\mbox{\boldmath$\bf x$})\\!\\!$ $\displaystyle=$
$\displaystyle\\!\\!{\mbox{\boldmath$\cal L$}}_{1}(\mbox{\boldmath$\bf
C$}_{1},\mbox{\boldmath$\bf u$}_{0})(\mbox{\boldmath$\bf x$}),$ (6.11)

which is equivalent to Eqs. (19.37) in [60]. Two terms in the right-hand side
of Eq. (6.111) can be expressed as ($v_{i}^{l}=v_{i}+v^{\Gamma-}_{i}$)

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!{\mbox{\boldmath$\cal
L$}}_{1}(\mbox{\boldmath$\bf C$}_{1},{\mbox{\boldmath$\bf
u$}}_{1})(\mbox{\boldmath$\bf x$})$ $\displaystyle=$
$\displaystyle{\mbox{\boldmath$\cal L$}}_{1}(\mbox{\boldmath$\bf
C$}^{0}_{1},V_{i}{\mbox{\boldmath$\bf u$}}_{1})(\mbox{\boldmath$\bf
x$})+\frac{1}{2}{\mbox{\boldmath$\cal L$}}_{1}(\mbox{\boldmath$\bf
C$}^{0}_{1},V_{i}^{\Gamma-}{\mbox{\boldmath$\bf u$}}_{1})(\mbox{\boldmath$\bf
x$}),$ (6.12)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!{\mbox{\boldmath$\cal
L$}}_{1}(\mbox{\boldmath$\bf C$}_{1},{\mbox{\boldmath$\bf
u$}}_{0})(\mbox{\boldmath$\bf x$})$ $\displaystyle=$
$\displaystyle{\mbox{\boldmath$\cal L$}}_{1}(\mbox{\boldmath$\bf
C$}^{0}_{1\epsilon},V_{i})(\mbox{\boldmath$\bf
x$})+\frac{1}{2}{\mbox{\boldmath$\cal L$}}_{1}(\mbox{\boldmath$\bf
C$}^{0}_{1\epsilon},V_{i}^{\Gamma-})(\mbox{\boldmath$\bf x$}).$ (6.13)

Thus, instead of the kernel $\mbox{\boldmath$\bf C$}(\mbox{\boldmath$\bf
x$},\mbox{\boldmath$\xi$})$ for an inhomogeneous medium (which is neither
convolution nor periodic in general), we introduced three convolution kernels
$\mbox{\boldmath$\bf C$}^{(0)}(\mbox{\boldmath$\xi$}),\mbox{\boldmath$\bf
C$}_{1}^{0}(\mbox{\boldmath$\xi$})$ and $\mbox{\boldmath$\bf
C$}_{1\epsilon}^{0}(\mbox{\boldmath$\xi$})$ (i.e., the obstacle O2) is
overcome) whereas other variables $\mbox{\boldmath$\bf
u$}_{1}(\mbox{\boldmath$\bf x$})$, $V_{i}{\mbox{\boldmath$\bf
u$}}_{1}(\mbox{\boldmath$\bf x$})$ and $V_{i}^{\Gamma-}{\mbox{\boldmath$\bf
u$}}_{1}(\mbox{\boldmath$\bf x$})$ are $\Omega_{00}$-periodic (i.e., the
obstacle O1) is overcome). Subsequent zero padding (6.3) (elimination of the
obstacle O3)) of both the kernels and integrands ensures the applicability of
the convolution theorem (6.32) that reduces the integral equation in the real
space (6.11) to the multiplication operations in the Fourier space

$\displaystyle\\!\\!\\!\widehat{\mbox{\boldmath$\bf
C$}}^{(0)}(\mbox{\boldmath$\bf k$})\widehat{\mbox{\boldmath$\bf
u$}}_{1}(\mbox{\boldmath$\bf k$})-\mbox{\boldmath$\bf
P$}^{(0)}\widehat{\mbox{\boldmath$\bf u$}}_{1}(\mbox{\boldmath$\bf
k$})=-\widehat{\mbox{\boldmath$\bf C$}}_{1}^{0}(\mbox{\boldmath$\bf
k$})\Big{[}\mbox{\boldmath$\cal F$}_{d}(V_{i}{\mbox{\boldmath$\bf
u$}})+\frac{1}{2}\mbox{\boldmath$\cal
F$}_{d}(V_{i}^{\Gamma-}{\mbox{\boldmath$\bf u$}})\Big{]}$ (6.14)
$\displaystyle+$ $\displaystyle\big{[}\mbox{\boldmath$\cal
F$}_{d}(\mbox{\boldmath$\bf P$}^{V}_{1}\mbox{\boldmath$\bf
u$}_{1})+\frac{1}{2}\mbox{\boldmath$\cal F$}_{d}(\mbox{\boldmath$\bf
P$}^{\Gamma-}_{1}\mbox{\boldmath$\bf
u$}_{1})\big{]}-\widehat{\mbox{\boldmath$\bf
C$}}_{1\epsilon}^{0}(\mbox{\boldmath$\bf k$})\big{[}\mbox{\boldmath$\cal
F$}_{d}(V_{i})+\frac{1}{2}\mbox{\boldmath$\cal
F$}_{d}(V_{i}^{\Gamma-})\big{]},$

where

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\mbox{\boldmath$\bf P$}^{(0)}$
$\displaystyle=$ $\displaystyle\int\mbox{\boldmath$\bf
C$}^{(0)}(\mbox{\boldmath$\bf x$},\hat{}\mbox{\boldmath$\bf
x$})~{}d\hat{}\mbox{\boldmath$\bf x$}=\mbox{\boldmath$\cal
F$}_{d}^{-1}\Big{(}\widehat{\mbox{\boldmath$\bf C$}}^{(0)}(\mbox{\boldmath$\bf
k$}){\bf 1}(\mbox{\boldmath$\bf k$})\Big{)},$ (6.15)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\mbox{\boldmath$\bf
P$}_{1}^{V}(\mbox{\boldmath$\bf x$})$ $\displaystyle=$
$\displaystyle\int\mbox{\boldmath$\bf C$}^{0}_{1}(\mbox{\boldmath$\bf
x$},\hat{}\mbox{\boldmath$\bf x$})V_{i}(\hat{}\mbox{\boldmath$\bf
x$})~{}d\hat{}\mbox{\boldmath$\bf x$}=\mbox{\boldmath$\cal
F$}_{d}^{-1}\Big{(}\widehat{\mbox{\boldmath$\bf
C$}}^{0}_{1}(\mbox{\boldmath$\bf k$})\widehat{V}_{i}(\mbox{\boldmath$\bf
k$})\Big{)},$ (6.16) $\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\mbox{\boldmath$\bf
P$}_{1}^{\Gamma-}(\mbox{\boldmath$\bf x$})$ $\displaystyle=$
$\displaystyle\int\mbox{\boldmath$\bf C$}^{0}_{1}(\mbox{\boldmath$\bf
x$},\hat{}\mbox{\boldmath$\bf x$})V^{\Gamma-}_{i}(\hat{}\mbox{\boldmath$\bf
x$})~{}d\hat{}\mbox{\boldmath$\bf x$}=\mbox{\boldmath$\cal
F$}_{d}^{-1}\Big{(}\widehat{\mbox{\boldmath$\bf
C$}}^{0}_{1}(\mbox{\boldmath$\bf
k$})\widehat{V}^{\Gamma-}_{i}(\mbox{\boldmath$\bf k$})\Big{)},$ (6.17)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!{\widehat{\mbox{\boldmath$\bf
P$}}_{1}}^{V}(\mbox{\boldmath$\bf k$})$ $\displaystyle=$
$\displaystyle\widehat{\mbox{\boldmath$\bf C$}}^{0}_{1}(\mbox{\boldmath$\bf
k$})\widehat{V}_{i}(\mbox{\boldmath$\bf k$}),\ \ \
{\widehat{\mbox{\boldmath$\bf P$}}_{1}}^{\Gamma-}(\mbox{\boldmath$\bf
k$})=\widehat{\mbox{\boldmath$\bf C$}}^{0}_{1}(\mbox{\boldmath$\bf
k$})\widehat{V}^{\Gamma-}_{i}(\mbox{\boldmath$\bf k$})$ (6.18)

and the indices N and $(\cdot)_{\rm e}$ are omitted for brevity in Eqs.
(6.14)-(6.18) as well as in the subsequent representations.

Applying iDFT to rearranged Eq. (6.14) leads to the implicit analog of the L-S
Eq. (6.7) (see for details [71], [72])

$\displaystyle\mbox{\boldmath$\bf u$}_{1}(\mbox{\boldmath$\bf x$})$
$\displaystyle=$ $\displaystyle\mbox{\boldmath$\bf g$}_{0}^{\bf
m}(\mbox{\boldmath$\bf x$})+{\mbox{\boldmath$\cal K$}}^{\bf
m}(\mbox{\boldmath$\bf u$}_{1})(\mbox{\boldmath$\bf x$})$ (6.19)

solved by the iteration method in accompanied with Eq. (6.14) (as in the pair
of Eqs. (6.7) and (6.8)). The convergence (algorithm stopping) criteria based
on either the compatibility condition or the imposed macroscopic condition
were considered [217] for locally elastic problems.

Comment 6.1. The scheme of the FFT method in both the LM and PM is presented
in Fig. 7. So, in preprocessing Blocks 1 and 2 (for the LM and PM,
respectively) we estimate the DFT of the following functions. According to the
properties of the DFT (6.3), in Block 1 (LM), we use either the analytical
representations of either integral operator (6.7) (see
$\widehat{}\mbox{\boldmath$\bf U$}(\mbox{\boldmath$\bf k$})$ in [221]) or the
differential equilibrium operator (2.4). In Block 2 (PM), the following
tensors are estimated one time (as in LM) in either the real space
($\mbox{\boldmath$\bf P$}_{1}^{V}(\mbox{\boldmath$\bf x$})$ and
$\mbox{\boldmath$\bf P$}_{1}^{\Gamma-}(\mbox{\boldmath$\bf x$})$) or frequency
domain ($\widehat{\mbox{\boldmath$\bf C$}}^{(0)}(\mbox{\boldmath$\bf
k$}),\widehat{\mbox{\boldmath$\bf C$}}_{1}^{(0)}(\mbox{\boldmath$\bf k$})$,
$\widehat{\mbox{\boldmath$\bf C$}}_{1\epsilon}^{(0)}(\mbox{\boldmath$\bf
k$})$, and $\widehat{V}_{i}(\mbox{\boldmath$\bf k$})$,
$\widehat{V}_{i}^{\Gamma-}(\mbox{\boldmath$\bf k$})$) which are also evaluated
by the FFT method. Blocks 3 and 4 are involved in the iteration scheme
corresponding to either Eqs. (6.7), (6.8) (in the LM) or (6.14), (6.19) (in
the PM).

Fig. 7: The scheme FFT methods in the LM and PM

Since Moulinec and Suquet [218] introduced the basic iteration FFT scheme for
Eqs. (6.7), (6.8), a family of improved accelerated solvers has been developed
to improve the convergence rate for highly-contrasted materials. These methods
(accelerated scheme [112], [208] and the augmented Lagrangian scheme [207])
and others proposed later [215], [217] share the same approach, which is based
on the introduction of a dual strain and/or dual stress-based formulation as
well as the modification of the polarization term in the L-S equation. It
would be interesting to generalize these approaches of the LM to their PM’s
counterparts.

All mentioned FFT schemes (based on the application of DFT (6.32) to L-S Eq.
(6.7)) have been presented for the strain field as the primary variable. An
alternative displacement-based FFT approach (based on the application of DFT
(6.31) to Eq. (2.4) (see [190], [191], [247]) involves solving directly for
the displacement field $\mbox{\boldmath$\bf u$}_{1}$, which replaces the
strain field as the primary variable. Then, rearrangement of the DFT form of
Eq. (2.4) and subsequent application of the iDFT leads to the operator
equation (which is not prescribed as a function of a matrix of coefficients)
in a real space

$\displaystyle\mbox{\boldmath$\cal A$}(\mbox{\boldmath$\bf
u$}_{1})=\mbox{\boldmath$\bf b$}^{\cal A},$ (6.20)

where the linear operator $\cal A$ defines the action on a vector
$\mbox{\boldmath$\bf u$}_{1}$. It is interesting that Eq. (6.19) in the PM is,
in fact, presented in the form (6.20), where

$\displaystyle\mbox{\boldmath$\cal A$}(\mbox{\boldmath$\bf
u$}_{1})=\mbox{\boldmath$\bf u$}_{1}(\mbox{\boldmath$\bf
x$})-\mbox{\boldmath$\cal K$}^{\bf\rm m}(\mbox{\boldmath$\bf
u$}_{1})(\mbox{\boldmath$\bf x$}),\ \ \ \mbox{\boldmath$\bf b$}^{\cal
A}(\mbox{\boldmath$\bf x$})={\mbox{\boldmath$\bf g$}}_{0}^{\bf\rm
m}(\mbox{\boldmath$\bf x$})$

and the displacement fluctuation $\mbox{\boldmath$\bf u$}_{1}$ is a primary
unknown variable in both the LM (2.4) and PM (6.11). Among the numerous
methods available for the solution of linear systems (6.20), Krylov method
(see the Krylov solver available in the PETSc library [28], see also [235])
such as the conjugate gradient (CG) method or the biconjugate gradient method
can be efficiently used. According [16] (see also [39]) all we need is to
provide the CG method by computing $\mbox{\boldmath$\cal
A$}(\mbox{\boldmath$\bf u$}_{1})$, where $\cal A$ is the linear operator
(rather than matrix). It was demonstrated [278] that the standard CG algorithm
in local micromechanics works well, even though the linear system is non-
symmetric (see also [171]). In the case of the LM, it is also demonstrated in
[278], [311] that contrary to the known accelerated FFT methods [112], [207]
improving the original basic scheme, the cost of one iteration remains
comparable to the basic scheme and the method proposed is insensitive to the
choice of the auxiliary reference medium ${}^{L}\\!\mbox{\boldmath$\bf
L$}^{(0)}$ (6.7).

Comment 6.2. The scheme presented in Fig. 7 can be applied to the different
methods (e.g. CG method) describing the interactions of Blocks 3 and 4 (rather
than only basic scheme [219]) for either LM or PM approaches. The specific
feature of the PM is preconditioning Block 2 analysing one inclusion (like
Block 1 Input in Fig. 2).

###  6.3  FFT of the nonlinear peridynamic equation

Equation (6.11) with the decomposition of the involved operators
${\mbox{\boldmath$\cal L$}}^{(0)}_{1}(\mbox{\boldmath$\bf
C$},{\mbox{\boldmath$\bf u$}})(\mbox{\boldmath$\bf x$})$ and
${\mbox{\boldmath$\cal L$}}_{1}(\mbox{\boldmath$\bf C$},{\mbox{\boldmath$\bf
u$}}_{0})(\mbox{\boldmath$\bf x$})$ can be presented in the following form.
So, the left-hand side of Eq. (6.11) is

$\displaystyle{\mbox{\boldmath$\cal L$}}^{(0)}(\mbox{\boldmath$\bf
C$}^{(0)},{\mbox{\boldmath$\bf u$}}_{1})(\mbox{\boldmath$\bf x$})$
$\displaystyle=$ $\displaystyle\int\mbox{\boldmath$\bf
C$}^{(0)}(\mbox{\boldmath$\bf x$}-\mbox{\boldmath$\bf
x$}^{\prime})(\mbox{\boldmath$\xi$}_{w\Gamma}+\mbox{\boldmath$\bf
u$}_{1}^{\prime}-\mbox{\boldmath$\bf u$}_{1})^{r}~{}dV_{\bf x^{\prime}},$
(6.21)

where $\mbox{\boldmath$\bf u$}_{1}=\mbox{\boldmath$\bf
u$}_{1}(\mbox{\boldmath$\bf x$})$, $\mbox{\boldmath$\bf
u$}_{1}^{\prime}=\mbox{\boldmath$\bf u$}_{1}(\mbox{\boldmath$\bf
x$}^{\prime})$,

$\mbox{\boldmath$\xi$}_{w\Gamma}=\mbox{\boldmath$\varepsilon$}^{w\Gamma}\mbox{\boldmath$\xi$}=\mbox{\boldmath$\varepsilon$}^{w\Gamma}(\mbox{\boldmath$\bf
x$}^{\prime}-\mbox{\boldmath$\bf x$}).$ (6.22)

By the use of the trinomial expansion

$\displaystyle(\mbox{\boldmath$\bf p$}+\mbox{\boldmath$\bf
q$}+\mbox{\boldmath$\bf
s$})^{r}=\sum_{a+b+c=r}\tbinom{n}{a,b,c}\mbox{\boldmath$\bf
p$}^{a}\mbox{\boldmath$\bf q$}^{b}\mbox{\boldmath$\bf s$}^{c},\ \ \
\tbinom{n}{a,b,c}:=\frac{n!}{a!\,b!\,c!}$ (6.23)

Eq. (6.21) is reduced to

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!{\mbox{\boldmath$\cal
L$}}^{(0)}(\mbox{\boldmath$\bf C$}^{(0)},{\mbox{\boldmath$\bf
u$}}_{1})(\mbox{\boldmath$\bf x$})$ $\displaystyle=$
$\displaystyle\\!\\!\\!\sum_{a+b+c=r}\\!\\!\\!\tbinom{r}{a,b,c}(-\mbox{\boldmath$\bf
u$}_{1})^{c}\\!\int\\!\\!\mbox{\boldmath$\bf C$}^{(0)}(\mbox{\boldmath$\bf
x$}-\mbox{\boldmath$\bf
x$}^{\prime})(\mbox{\boldmath$\xi$}_{w\Gamma})^{a}(\mbox{\boldmath$\bf
u$}_{1}^{\prime})^{b}~{}dV_{\bf x^{\prime}}.$ (6.24)

Similarly,

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!{\mbox{\boldmath$\cal
L$}}_{1}(\mbox{\boldmath$\bf C$}^{0}_{1},{\mbox{\boldmath$\bf
u$}}_{1})(\mbox{\boldmath$\bf x$})$ $\displaystyle=\\!\\!\\!\\!\\!\\!$
$\displaystyle\sum_{a+b+c=r}\\!\\!\\!\\!\\!\\!\tbinom{r}{a,b,c}(-\mbox{\boldmath$\bf
u$}_{1})^{c}\\!\\!\int\\!\mbox{\boldmath$\bf C$}^{0}_{1}(\mbox{\boldmath$\bf
x$}-\mbox{\boldmath$\bf
x$}^{\prime})(\mbox{\boldmath$\xi$}_{w\Gamma})^{a}(\mbox{\boldmath$\bf
u$}_{1}^{\prime})^{b}V_{i}^{l\Gamma}\\!(\mbox{\boldmath$\bf
x$}^{\prime})dV_{\bf x^{\prime}}.$ (6.25)

Solution of PD Eq. (6.21) with polynomial nonlinearity of the integrand was
analyzed in [188] and [189] (see also [84]) by FFT method for the homogeneous
matrix $\mbox{\boldmath$\bf C$}_{1}^{0}\equiv{\bf 0}$. We use a similar
decomposition (6.25) and (6.3) for composite material ($\mbox{\boldmath$\bf
C$}_{1}^{0}\not\equiv{\bf 0}$).

We introduce the new convolution kernels

$\displaystyle\mbox{\boldmath$\bf C$}^{(0)\xi a}(\mbox{\boldmath$\bf
x$}-\mbox{\boldmath$\bf x$}^{\prime})$ $\displaystyle=$
$\displaystyle\mbox{\boldmath$\bf C$}^{(0)}(\mbox{\boldmath$\bf
x$}-\mbox{\boldmath$\bf x$}^{\prime})(\mbox{\boldmath$\xi$}_{w\Gamma})^{a},$
$\displaystyle\mbox{\boldmath$\bf C$}^{0\xi a}_{1}(\mbox{\boldmath$\bf
x$}-\mbox{\boldmath$\bf x$}^{\prime})$ $\displaystyle=$
$\displaystyle\mbox{\boldmath$\bf C$}^{0}_{1}(\mbox{\boldmath$\bf
x$}-\mbox{\boldmath$\bf x$}^{\prime})(\mbox{\boldmath$\xi$}_{w\Gamma})^{a},$
(6.26)

with initial values

$\displaystyle\mbox{\boldmath$\bf C$}^{(0)\xi 0}(\mbox{\boldmath$\bf
x$}-\mbox{\boldmath$\bf x$}^{\prime})=\mbox{\boldmath$\bf
C$}^{(0)}(\mbox{\boldmath$\bf x$}-\mbox{\boldmath$\bf x$}^{\prime}),\ \ \
\mbox{\boldmath$\bf C$}^{0\xi 0}_{1}(\mbox{\boldmath$\bf
x$}-\mbox{\boldmath$\bf x$}^{\prime})=\mbox{\boldmath$\bf
C$}^{0}_{1}(\mbox{\boldmath$\bf x$}-\mbox{\boldmath$\bf x$}^{\prime})$ (6.27)

that simplifies Eqs. (6.24) and (6.25)

$\displaystyle\sum_{a+b+c=r}\tbinom{r}{a,b,c}(-\mbox{\boldmath$\bf
u$}_{1})^{c}\int\mbox{\boldmath$\bf C$}^{(0)\xi a}(\mbox{\boldmath$\bf
x$}-\mbox{\boldmath$\bf x$}^{\prime})(\mbox{\boldmath$\bf
u$}_{1}^{\prime})^{b}~{}dV_{\bf x^{\prime}}$ (6.28) $\displaystyle=$
$\displaystyle-\sum_{a+b+c=r}\tbinom{r}{a,b,c}(-\mbox{\boldmath$\bf
u$}_{1})^{c}\int\mbox{\boldmath$\bf C$}^{0\xi a}_{1}(\mbox{\boldmath$\bf
x$}-\mbox{\boldmath$\bf x$}^{\prime})(\mbox{\boldmath$\bf
u$}_{1}^{\prime})^{b}V_{i}^{l\Gamma}(\mbox{\boldmath$\bf
x$}^{\prime})~{}dV_{\bf x^{\prime}}$

Zero-pudding of both the convolution kernels and prime variable in Eq. (6.28)
transforms the convolution operations into circular convolutions:

$\displaystyle\\!\\!\\!\sum_{a+b+c=r}\tbinom{r}{a,b,c}(-\mbox{\boldmath$\bf
u$}_{1})^{c}\mbox{\boldmath$\bf C$}^{(0)\xi a}\odot(\mbox{\boldmath$\bf
u$}_{1}^{\prime})^{b}(\mbox{\boldmath$\bf x$})$ (6.29) $\displaystyle=$
$\displaystyle\sum_{a+b+c=r}\tbinom{r}{a,b,c}(-\mbox{\boldmath$\bf
u$}_{1})^{c}\mbox{\boldmath$\bf C$}^{0\xi
a}_{1}\odot\big{[}(\mbox{\boldmath$\bf
u$}_{1}^{\prime})^{b}V_{i}^{l\Gamma}\big{]}(\mbox{\boldmath$\bf x$}),$

hereafter the indices N and $(\cdot)_{\rm e}$ (see the notations in Section 4)
are omitted for brevity. Thanks to the Convolution Theorem (4.10), we can
compute the circular convolution in (6.28) by means of the inverse Fourier
transform $\mbox{$\cal F$}_{d}^{-1}$ of the product of Fourier transforms:

$\displaystyle\mbox{\boldmath$\bf C$}^{(0)\xi a}\odot\mbox{\boldmath$\bf
u$}_{1}^{b}(\mbox{\boldmath$\bf x$})$ $\displaystyle=$
$\displaystyle\mbox{$\cal F$}_{d}^{-1}\big{[}\mbox{$\cal
F$}_{d}(\mbox{\boldmath$\bf C$}^{(0)\xi a})(\mbox{\boldmath$\bf
k$})\mbox{$\cal F$}_{d}(\mbox{\boldmath$\bf u$}_{1}^{b})(\mbox{\boldmath$\bf
k$})\big{]},$ (6.30) $\displaystyle\mbox{\boldmath$\bf C$}_{1}^{0\xi
a}\odot(V_{i}^{l\Gamma}\mbox{\boldmath$\bf u$}_{1}^{b})(\mbox{\boldmath$\bf
x$})$ $\displaystyle=$ $\displaystyle\mbox{$\cal
F$}_{d}^{-1}\big{[}\mbox{$\cal F$}_{d}(\mbox{\boldmath$\bf C$}_{1}^{0\xi
a})(\mbox{\boldmath$\bf k$})\mbox{$\cal
F$}_{d}(V_{i}^{l\Gamma}\mbox{\boldmath$\bf u$}_{1}^{b})(\mbox{\boldmath$\bf
k$})\big{]}.$ (6.31)

Additionally, according to the Inverse Theorem, we obtain

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!(\mbox{\boldmath$\bf
u$}_{1})^{c}\big{(}\mbox{\boldmath$\bf C$}^{(0)\xi a}\odot\mbox{\boldmath$\bf
u$}_{1}^{b}\big{)}(\mbox{\boldmath$\bf x$})$ $\displaystyle=$
$\displaystyle\mbox{$\cal F$}_{d}^{-1}\Big{\\{}\mbox{$\cal
F$}_{d}(\mbox{\boldmath$\bf u$}_{1}^{c})(\mbox{\boldmath$\bf k$})$ (6.32)
$\displaystyle\odot$ $\displaystyle\big{[}\mbox{$\cal
F$}_{d}(\mbox{\boldmath$\bf C$}^{(0)\xi a})(\mbox{\boldmath$\bf
k$})\mbox{$\cal F$}_{d}((\mbox{\boldmath$\bf
u$}_{1}^{\prime})^{b})(\mbox{\boldmath$\bf k$})\big{]}\Big{\\}},$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!(\mbox{\boldmath$\bf
u$}_{1})^{c}\big{(}\mbox{\boldmath$\bf C$}_{1}^{0\xi
a}\odot(V_{i}^{l\Gamma}(\mbox{\boldmath$\bf
u$}_{1}^{\prime})^{b}\big{)}(\mbox{\boldmath$\bf x$})$ $\displaystyle=$
$\displaystyle\mbox{$\cal F$}_{d}^{-1}\Big{\\{}\mbox{$\cal
F$}_{d}(\mbox{\boldmath$\bf u$}_{1}^{c})(\mbox{\boldmath$\bf k$})$ (6.33)
$\displaystyle\odot$ $\displaystyle\big{[}\mbox{$\cal
F$}_{d}(\mbox{\boldmath$\bf C$}_{1}^{0\xi a})(\mbox{\boldmath$\bf
k$})\mbox{$\cal F$}_{d}(V_{i}^{l\Gamma}(\mbox{\boldmath$\bf
u$}_{1}^{\prime})^{b})(\mbox{\boldmath$\bf k$})\big{]}\Big{\\}}.$

In particular, by the DFT/iDFT methods, we can estimate the terms of Eq.
(6.28) corresponding to $b=0$

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\mbox{\boldmath$\bf P$}^{(0)\xi
a}(\mbox{\boldmath$\bf x$})$ $\displaystyle:=$
$\displaystyle\int\mbox{\boldmath$\bf C$}^{(0)\xi a}(\mbox{\boldmath$\bf
x$}-\mbox{\boldmath$\bf x$}^{\prime})~{}d\mbox{\boldmath$\bf
x$}^{\prime}=\mbox{$\cal F$}^{-1}\Big{(}\widehat{\mbox{\boldmath$\bf
C$}}^{(0)\xi a}(\mbox{\boldmath$\bf k$})\widehat{\bf 1}(\mbox{\boldmath$\bf
k$})\Big{)},$ (6.34)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\mbox{\boldmath$\bf
P$}_{1}^{V\xi a}(\mbox{\boldmath$\bf x$})$ $\displaystyle:=$
$\displaystyle\int\mbox{\boldmath$\bf C$}_{1}^{0\xi a}(\mbox{\boldmath$\bf
x$}-\mbox{\boldmath$\bf x$}^{\prime})V_{i}^{l\Gamma}(\mbox{\boldmath$\bf
x$}^{\prime})~{}d\mbox{\boldmath$\bf x$}^{\prime}=\mbox{$\cal
F$}^{-1}\Big{(}\widehat{\mbox{\boldmath$\bf C$}}_{1}^{0\xi
a}(\mbox{\boldmath$\bf k$})\widehat{V}_{i}^{l\Gamma}(\mbox{\boldmath$\bf
k$})\Big{)},$ (6.35)

which the generalization of the terms (6.15) and (6.16), respectively.
Substitution of Eqs. (6.34)-(6.35) into Eq. (6.29) leads to the representation
generalizing Eq. (6.19). However, a more detailed analysis of this approach is
beyond the scope of the current study and deserve further investigations.

Comment 6.3. We considered in Section 7 only one type of nonlinear
peridynamics, precisely nonlinear elasticity (6.21). In the LM over the past
30 years (after the seminal paper [218]), several high-performance FFT
algorithms with improved convergence rates and solution accuracy have been
developed to study the non-linear response (such as plasticity,
viscoplasticity, damage, fracture, and fatigue, etc.) of microstructures with
arbitrary phase contrast subjected to finite deformations (see the
comprehensive reviews [83], [129], [191], [247], [248]), which can be computed
in desktop computers without the need to use high-performance computational
clusters or supercomputers. They have also been exploited to study multi-
physics systems involving the coupled mechanical, thermal, electrical,
magnetic, and pyroelectric fields; for example, the chemo-thermo-mechanical
study of batteries, including phase-field damage, which may rely on different
fields coupled by their corresponding PDEs. All mentioned methods and problems
are looking forward to their generalizing to the PM.

##  7 Reduced order models in peridynamic micromechanicss

###  7.1 Clustering discretization method in local micromechanics

In local micromechanics, for the periodic structure CMs, the total fields
(displacement, strain, or stress) caused in $v^{(s)}$ by the uniform loading
fields either
$\mbox{\boldmath$\varepsilon$}^{w_{\Gamma}}=\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}=$const.
or
$\mbox{\boldmath$\sigma$}^{w_{\Gamma}}=\mbox{$\langle$}\mbox{\boldmath$\sigma$}\mbox{$\rangle$}=$const.
(2.3) (at $l_{\delta}/a\to 0$) and the transformation fields (2.2) are
represented in the following form $(\mbox{\boldmath$\bf x$}\in v^{(s)}$)

$\displaystyle\mbox{\boldmath$\vartheta$}(\mbox{\boldmath$\bf
z$})=\mbox{\boldmath$\bf A$}^{(s)}_{\theta\zeta}(\mbox{\boldmath$\bf
z$})\mbox{$\langle$}\mbox{\boldmath$\zeta$}\mbox{$\rangle$}+\sum_{r=0}^{N}\mbox{\boldmath$\bf
D$}^{sr}_{\theta}(\mbox{\boldmath$\bf z$}){\mbox{\boldmath$\beta$}}^{(r)},$
(7.1)

which is a slight modification of [105]. For compactness of representations
(7.1), one introduces the substitutions

$(\mbox{\boldmath$\bf
u$}_{1},\mbox{\boldmath$\sigma$},\mbox{\boldmath$\eta$})\leftrightarrow\mbox{\boldmath$\vartheta$},\
\ \
(\mbox{\boldmath$\varepsilon$},\mbox{\boldmath$\sigma$})\leftrightarrow\mbox{\boldmath$\zeta$},\
\ \ [\mbox{\boldmath$\bf x$},(\hat{}\mbox{\boldmath$\bf
x$},\mbox{\boldmath$\bf x$})]\leftrightarrow\mbox{\boldmath$\bf z$}.$ (7.2)

Here, the $\mbox{\boldmath$\bf A$}_{\theta\zeta}(\mbox{\boldmath$\bf z$})$ is
the mechanical influence function (do not miss with the influence functions in
peridynamics, see [252]) whereas he $\mbox{\boldmath$\bf
D$}^{sr}_{\theta}(\mbox{\boldmath$\bf z$})$ is the eigenfield influence
functions.

In the offline stage, the high-fidelity UC is represented by $n$ voxels or
elements denoted as $\omega_{m}$ ($m=1,\ldots,n)$. K-means clustering gropes
data vectors $\mbox{\boldmath$\bf A$}_{\sigma\varepsilon m}$ into a predefined
number $K$ of clusters $\Omega=\Omega^{1}\cup\ldots\cup\Omega^{K}$, based on
Euclidean distance as a similarity measure; any points $\mbox{\boldmath$\bf
x$}\in\Omega^{J}$ ($J=1,\ldots,K$) have the same response
$\mbox{\boldmath$\sigma$}(\mbox{\boldmath$\bf x$})$ close to each other. Owing
to a large number of data points, an unsupervised machine learning technique
such as $k$-means or Self-Organizing Map, SOM, (see [153], [184]) is adopted
to perform the decomposition of UC. Combined with the stress concentration
tensor, the mathematic definition of $k$‐means clustering is based on
optimization of a criterium function for the $k$ sets $\mbox{\boldmath$\bf
S$}=(\Omega^{1},\ldots,\Omega^{K})$

$\mbox{\boldmath$\bf S$}={\rm argmin}_{{}_{\mbox{\boldmath$\bf
S$}^{\prime}}}~{}~{}~{}{\rm\bf SSE},\ \ \ {\rm\bf
SSE}:=\sum_{J=1}^{K}\sum_{m\in S^{J}}||\mbox{\boldmath$\bf
A$}_{\sigma\varepsilon m}-\mbox{\boldmath$\bf
A$}_{\sigma\varepsilon}^{J}||^{2},$ (7.3)

where $\mbox{\boldmath$\bf A$}_{\sigma\varepsilon}^{J}$ is the average of
strain concentration factor in the $J$th cluster $\Omega^{J}$
$\mbox{$\langle$}\mbox{\boldmath$\sigma$}\mbox{$\rangle$}_{J}=\mbox{\boldmath$\bf
A$}_{\sigma\varepsilon}^{I}\mbox{\boldmath$\varepsilon$}^{w\Gamma}$ whereas
the norm $||(\cdot)||$ is defined as e.g., Frobenius norm for matrix.

Averaging of Eq. (7.1) over the domain $\Omega^{J}$ yields

$\displaystyle\mbox{\boldmath$\sigma$}_{I}$ $\displaystyle=$
$\displaystyle\mbox{\boldmath$\bf
A$}_{\sigma\varepsilon}^{I}\mbox{\boldmath$\varepsilon$}^{w\Gamma}+\sum_{J=1}^{K}\mbox{\boldmath$\bf
D$}^{IJ}_{\sigma}\mbox{\boldmath$\beta$}^{J},$ (7.4)

where the physical meaning of ${\mbox{\boldmath$\bf D$}}^{JL}_{\sigma}$ is the
local stress field $\mbox{\boldmath$\sigma$}^{J}$ in the domain
$\mbox{\boldmath$\bf x$}\in\Omega^{J}$ at
$\mbox{$\langle$}\mbox{\boldmath$\varepsilon$}\mbox{$\rangle$}=\mbox{\boldmath$\varepsilon$}^{w\Gamma}\equiv{\bf
0}$ and the uniform eigenstrain of unit magnitude prescribed in every single
domain $\Omega^{J}$, while all other domains are free of any transformation.

In nonlinear micromechanics, the strains are assumed to be small, but both
elastic and inelastic behavior is that at each instant of loading the additive
decomposition takes place

$\displaystyle\mbox{\boldmath$\sigma$}_{I}(\mbox{\boldmath$\bf x$})$
$\displaystyle=$ $\displaystyle\mbox{\boldmath$\sigma$}_{I}^{\rm
lin}(\mbox{\boldmath$\bf x$})+\mbox{\boldmath$\sigma$}_{I}^{\rm
re}(\mbox{\boldmath$\bf x$}),\ \ {\rm or}\ \ \
\mbox{\boldmath$\varepsilon$}_{I}(\mbox{\boldmath$\bf
x$})=\mbox{\boldmath$\varepsilon$}_{I}^{\rm lin}(\mbox{\boldmath$\bf
x$})+\mbox{\boldmath$\varepsilon$}_{I}^{\rm in}(\mbox{\boldmath$\bf x$}).$
(7.5)

Here $\mbox{\boldmath$\sigma$}_{I}^{\rm lin}(\mbox{\boldmath$\bf x$})$ and
$\mbox{\boldmath$\varepsilon$}_{I}^{\rm lin}(\mbox{\boldmath$\bf x$})$ are the
linear elastic fields corresponding to the constitutive relations (2.2) and
the ineelastic fields (analog of the eigenstress
$\mbox{\boldmath$\alpha$}(\mbox{\boldmath$\bf x$})$ and eigenstrain
$\mbox{\boldmath$\beta$}(\mbox{\boldmath$\bf x$})$) are linked as

$\displaystyle\mbox{\boldmath$\varepsilon$}_{I}^{\rm in}(\mbox{\boldmath$\bf
x$})$ $\displaystyle=$ $\displaystyle-\mbox{\boldmath$\bf
M$}^{I}\mbox{\boldmath$\sigma$}_{I}^{\rm re}(\mbox{\boldmath$\bf x$}),\ \
\mbox{\boldmath$\sigma$}_{I}^{\rm re}(\mbox{\boldmath$\bf
x$})=-\mbox{\boldmath$\bf L$}^{I}\mbox{\boldmath$\varepsilon$}_{I}^{\rm
in}(\mbox{\boldmath$\bf x$}).$ (7.6)

We assume that clustering is performed and presented for piecewise constant
fields $\mbox{\boldmath$\varepsilon$}_{I},\mbox{\boldmath$\sigma$}_{I}$ and
$\mbox{\boldmath$\varepsilon$}_{J}^{\rm in}$, where
$\mbox{\boldmath$\varepsilon$}_{J}^{\rm in}$ is substituted instead of
eigenstrain $\mbox{\boldmath$\beta$}^{J}$. Many inelastic constitutive laws
relate either $\mbox{\boldmath$\sigma$}_{J}^{\rm re}$ to
$\mbox{\boldmath$\varepsilon$}_{J}$ or $\mbox{\boldmath$\varepsilon$}_{J}^{\rm
in}$ to $\mbox{\boldmath$\sigma$}_{J}$: $\mbox{\boldmath$\sigma$}_{J}^{\rm
re}(\mbox{\boldmath$\bf
x$})={\mbox{\boldmath$\mathfrak{g}$}}(\mbox{\boldmath$\varepsilon$}_{J}),\ \ \
\mbox{\boldmath$\varepsilon$}_{J}^{\rm
in}={\mbox{\boldmath$\mathfrak{f}$}}(\mbox{\boldmath$\sigma$}_{J})$. It
provides the following systems of governing equations for the evaluation of
local fields $(I=1,\ldots,K$)

$\displaystyle\mbox{\boldmath$\varepsilon$}_{I}=\mbox{\boldmath$\bf
A$}_{\varepsilon\varepsilon}^{I}\mbox{\boldmath$\varepsilon$}^{w\Gamma}-\sum_{J=1}^{K}\mbox{\boldmath$\bf
D$}^{IJ}_{\varepsilon}\mbox{\boldmath$\bf
M$}^{J}{\mbox{\boldmath$\mathfrak{g}$}}(\mbox{\boldmath$\varepsilon$}_{J}),\ \
\ \mbox{\boldmath$\sigma$}_{I}=\mbox{\boldmath$\bf
A$}_{\sigma\varepsilon}^{I}\mbox{\boldmath$\varepsilon$}^{w\Gamma}-\sum_{J=1}^{K}\mbox{\boldmath$\bf
D$}^{IJ}_{\sigma}\mbox{\boldmath$\mathfrak{f}$}(\mbox{\boldmath$\sigma$}_{J}).$
(7.7)

Numerical solution of Eqs. (7.7) (and their incremental analogs) are usually
found by Newton’s iteration method (see, e.g., [184]). In so doing, analysis
of deformation of each cluster $\Omega^{J}$ ($J=1,\ldots,M$) (possibly