2.5cm2.5cm2cm2cm

# Deep material networks for fiber suspensions with infinite material contrast

Benedikt Sterr  Karlsruhe Institute of Technology (KIT), Institute of
Engineering Mechanics  Sebastian Gajek  Karlsruhe Institute of Technology
(KIT), Institute of Engineering Mechanics  Andrew Hrymak  Matti Schneider
Thomas Böhlke  Karlsruhe Institute of Technology (KIT), Institute of
Engineering Mechanics

###### Abstract

We extend the laminate based framework of direct Deep Material Networks (DMNs)
to treat suspensions of rigid fibers in a non-Newtonian solvent. To do so, we
derive two-phase homogenization blocks that are capable of treating
incompressible fluid phases and infinite material contrast. In particular, we
leverage existing results for linear elastic laminates to identify closed form
expressions for the linear homogenization functions of two-phase layered
emulsions. To treat infinite material contrast, we rely on the repeated
layering of two-phase layered emulsions in the form of coated layered
materials. We derive necessary and sufficient conditions which ensure that the
effective properties of coated layered materials with incompressible phases
are non-singular, even if one of the phases is rigid.
With the derived homogenization blocks and non-singularity conditions at hand,
we present a novel DMN architecture, which we name the Flexible DMN (FDMN)
architecture. We build and train FDMNs to predict the effective stress
response of shear-thinning fiber suspensions with a Cross-type matrix
material. For 31 fiber orientation states, six load cases, and over a wide
range of shear rates relevant to engineering processes, the FDMNs achieve
validation errors below 4.31% when compared to direct numerical simulations
with Fast-Fourier-Transform based computational techniques. Compared to a
conventional machine learning approach introduced previously by the consortium
of authors, FDMNs offer better accuracy at an increased computational cost for
the considered material and flow scenarios.

Keywords: Deep Material Network; Fiber-reinforced composites; Non-Newtonian
suspension; Infinite material contrast; Supervised machine learning; Effective
viscosity

##  1 Introduction

###  1.1 State of the art

Especially because of their high specific stiffness, fiber reinforced
composites are widely used for lightweight design, particularly in the
transport and energy sectors [1, 2]. During the design process of fiber
reinforced composite components, combinations of molding simulations and
structural simulations are central to facilitate use case appropriate and
material oriented design through digital twins [3] and virtual process chains
[4, 5, 6]. In molding simulations of fiber reinforced composite parts,
accurate prediction of the local fiber suspension viscosity is vital to
correctly estimate process parameters [7, 8] and other engineering quantities.
In particular, the fiber orientation and the fiber volume distributions [9],
the flow fields, [10], as well as the final part properties [11, 12, 10]
depend on the suspension viscosity. However, in component scale simulations,
it is computationally infeasible to fully resolve the microstructure of the
whole domain of interest to compute the suspension viscosity. For this
purpose, analytical and computational homogenization methods are valuable
tools to provide viscosity estimates for molding simulations. However,
significant challenges arise for the holistic analytical modelling of the
suspension viscosity, because the suspension viscosity depends on the local
microstructure [13], the fiber volume fraction [14], and the fiber orientation
state [15]. Additionally, the flow field [16, 17] and the melt temperature
[18] influence the suspension viscosity as well. Furthermore, typical matrix
materials in fiber suspensions show non-Newtonian behavior, which adds
additional complexity to the derivation of appropriate models for the
suspension viscosity.

Depending on the fiber concentration of the suspension and the behavior of the
matrix material, different homogenization approaches have been proposed.
According to Pipes et al. [19, §1] the fiber concentration regimes may be
defined using the fiber volume fraction $c_{\sf{F}}$ and the fiber aspect
ratio $r_{\sf{a}}$ as follows. The dilute regime is defined through
${c_{\sf{F}}<(1/r_{\sf{a}})^{2}}$, the semi-dilute regime through
${(1/r_{\sf{a}})^{2}<c_{\sf{F}}<1/r_{\sf{a}}}$, the concentrated regime
through ${1/r_{\sf{a}}<c_{\sf{F}}}$, and the hyperconcentrated regime through
${r_{\sf{a}}>100}$ [19]. Goddard [20, 21], and Souloumiac and Vincent [22]
proposed self-consistent analytical models for different concentration
regimes. The models by Goddard are restricted to the dilute and semi-dilute
regimes, while the models by Souloumiac and Vincent are applicable to the
dilute, semi-dilute, and concentrated regimes. Both models agree well with
experimental results [23] qualitatively. However, the quantitative prediction
accuracy could be improved and varies strongly depending on the shear rate and
fiber volume fraction. Similarly, a semi-analytical model proposed by Férec et
al. [24] for Ellis and Carreau-type matrix behavior replicates the steady
state flow solution of a simple shear flow, but its accuracy varies with the
shear rate and the fiber volume fraction. By incorporating uniformly
distributed fiber misalignments with an orientation averaging approach [25,
26], Pipes et al. [19, 27, 28, 29] derived models for the concentrated and
hyper-concentrated regime. The models are only applicable to collimated fiber
arrays, but their predictions agree well with experimental studies by Binding
[30]. Combining orientation averaging with a deformation mode and
microstructure dependent informed isotropic viscosity, Favaloro et al. [31]
derived a model that aims to improve available molding simulation solvers with
only small modifications. Depending on the deformation mode and the
approximated anisotropic viscosity, the prediction quality of the model
varies. However, molding simulations using the model successfully predicted
the shell-core effect that is common in a wide variety of fiber molding
applications [9, 11]. Overall, predicting the suspension viscosity accurately
over the wide variety of processing conditions used for engineering systems
proves difficult, which includes recent developments [32, 33].

The development of homogenization methods is accompanied by considerable
difficulties in studying the fiber suspension viscosity experimentally. In
rheometer studies, certain transient effects like fiber breakage and fiber
reorientation are difficult to quantify during the experiment. Thus, it is
hard to link a specific fiber orientation state and a specific fiber length
distribution to a measured viscosity [30, 34]. Furthermore, the measurement
results may be affected by the interaction between the fibers and the
measurement devices [35].

In view of these difficulties, computational techniques offer an extraordinary
approach to study the effective viscosity of fiber suspensions with non-
Newtonian solvents by resolving the local fields of interest and providing
insights, which are hard to obtain otherwise. Combining a mass tracking
algorithm for the free surface representation, a lattice Boltzmann method for
fluid flow, and an immersed boundary procedure for the interaction between
fluid and rigid particles, Švec et al. [36] studied slump tests of rigid
fibers and rigid spherical particles suspended in a Bingham-type fluid. They
compared slump tests of the suspension with a slump test of the pure matrix
material, and observed a smaller spread and an increased height in the test of
the suspension material, which implies an increased effective yield stress in
the suspension. Domurath et al. [37] employed a Finite Element Method (FEM)
based approach to study the transversely isotropic fluid equation by Ericksen
[25]. They investigated the rheological coefficients and found that the model
by Souloumiac and Vincent [22] overpredicts the orientation dependence of a
rheological coefficient. Extending work by Bertóti et al. [38] on suspensions
with Newtonian solvents, Sterr et al. [17] used Fast Fourier Transform (FFT)
based computational techniques [39] and the RVE method [40] to study the
effective viscosity of fiber suspensions with non-Newtonian solvents. They
investigated the effects of anisotropic shear-thinning on the effective
suspension viscosity for varying fiber volume fractions, shear rates, and flow
scenarios.

Accurate computational predictions of the microscopic behavior of a composite
often require significant computational effort [41, 42, 43, 44, 45], which
motivates the use of fast data based surrogate models for component scale
simulations [46]. Ashwin et al. [47, 48] trained a Multi-Layer-Perceptron, a
Convolutional Neural Network, and a U-Net [49] on particle resolved simulation
data to predict fluid forces in dense ellipsoidal particle suspensions. They
restricted to Newtonian matrix behavior and trained the networks on data for
various Reynolds numbers and fiber volume fractions. Depending on the
considered combination of Reynolds number and particle volume fraction, the
prediction accuracy of the surrogate models varies. Boodaghidizaji et al. [50]
use a multi-fidelity approach with neural networks and Gaussian processes to
predict the steady state viscosity of fiber suspensions with a Newtonian
solvent. To form the training data set, they combine low-fidelity estimates
from constitutive equations with high-fidelity data obtained from numerically
solving the involved partial differential equation system. The prediction
accuracy of the multi-fidelity neural network and the multi-fidelity Gaussian
process for simple shear flow depends strongly on the investigated parameters,
especially the fiber volume fraction. Sterr et al. [51] derived four models
for fiber suspensions with a Cross-type matrix fluid by combining FFT-based
computational homogenization techniques with supervised machine learning. They
investigated the anisotropic shear-thinning characteristics of the suspension
viscosity for a variety of fiber orientation states via computational
experiments, and formulated model candidates based on the observed phenomena.
Using supervised machine learning techniques, they identified the model
parameters from computational data, so that three of the four models were able
to predict the fiber suspension viscosities to engineering accuracy for a wide
range of engineering load cases. Generally, in addition to accurate
estimation, extrapolation beyond the training data, as well as ensuring
thermodynamical consistency prove challenging in the construction of surrogate
models.

In the context of solid materials without kinematic constraints, Liu et
al.[52, 53] proposed Deep Material Networks (DMNs) as surrogate models for the
full-field computational homogenization of microstructured materials. Their
approach is based on nesting rotated laminates in an $N$-ary tree structure of
$N$-phase laminates, and thus constructing a micromechanically motivated deep
learning architecture. The volume fractions and rotations of the DMN are then
identified via supervised machine learning on linear elastic data. Remarkably,
even if a DMN is trained on data obtained by solving linear homogenization
problems, the predictions of the DMN for non-linear homogenization problems
are impressively accurate. Gajek et al. [54] further developed the DMN
architecture into the rotation free direct DMN architecture. Direct DMNs
feature a faster and more robust training process, as well as an efficient
evaluation scheme for non-linear problems. Gajek et al. [54] also showed that
the thermodynamic consistency of the laminates is preserved in (direct) DMNs,
such that the resulting DMN is thermodynamically consistent as well.
Furthermore, Gajek et al. [54] proved that non-linear homogenization is
determined by linear homogenization to first order for two-phase materials. In
a later article, Gajek et al. [46] used direct DMNs to accelerate two-scale FE
simulations of fiber reinforced composites by augmenting direct DMNs with the
fiber orientation interpolation concept introduced by Köbler et al. [55].
Alternatively, DMN parameters may be interpolated with regard to
microstructural parameters by a posteriori interpolation as proposed by Liu et
al. [56] and Huang et al. [57], or by augmenting DMNs with neural networks as
proposed by Li [58]. To leverage the capabilities of direct DMNs in concurrent
thermomechanical two-scale simulations of composite components, Gajek et al.
[59] further extended the direct DMN architecture to incorporate
thermomechanical coupling. Also in a thermomechanical setting, Shin et al.
[60] trained DMNs on linear thermoelastic data instead of linear elastic data,
which improved the quality of fit for the effective thermal expansion
properties, but only slightly affected the non-linear prediction quality of
the DMNs. Additionally, they employed DMNs for uncertainty quantification, and
for the inverse problem of optimizing a thermal boundary condition to achieve
a desired thermo-elasto-viscoplastic response. By developing an inelastically-
informed training strategy for DMNs, Dey et al. [61] successfully predicted
the creep behavior of fiber reinforced thermoplastics, which involves multiple
scales in both space and time. This enabled the inverse calibration of
parameters for creep and plasticity constitutive equations by using DMNs as
surrogates for otherwise costly FFT-based computations [62]. Furthermore, Dey
et al. [63] leveraged DMNs in combination with fiber orientation interpolation
[55] to characterize the behavior of fiber reinforced thermoplastics including
damage, plasticity, and creep. Overall, DMNs were extended and applied to
treat a wide variety of problems, such as interface damage [64], the modeling
of multiscale strain localization [65], problems involving woven materials
[66] and porous materials [67], as well as the architecture independent
treatment of multi-phase composites [68].

For composites with kinematically unconstrained solid phases, the presented
literature provides compelling evidence for the accuracy and versatility of
(direct) DMNs. However, the (direct) DMN architecture has yet to be applied to
composites involving fluids or infinite material contrast, which is required
for suspensions of rigid fibers. Because (direct) DMNs are based on laminates
of solids, the treatment of fluids with a DMN architecture requires a
different type of building block. Additionally, treatment of infinite material
contrast with DMNs is challenging, because the rank-one laminates of a DMN
have singular effective properties if one of the phases is singular, i.e.,
rigid. Therefore, the singular effective properties of the rank-one laminates
may propagate through the whole DMN during the training and evaluation
processes, rendering the effective properties of the DMN singular as well. In
this article we address these issues and propose an architecture to treat
suspensions of rigid fibers.

###  1.2 Contributions

In this work, we extend the direct DMN architecture for kinematically
unconstrained solid phases with finite material contrast to an architecture
that is able to treat suspensions of rigid fibers. In particular, this task
requires the treatment of fluid phases and infinite material contrast. We name
this architecture the Flexible DMN (FDMN) architecture. To treat fluid
materials with FDMNs, we derive homogenization blocks for layered emulsions
that are governed by Stokes flow, and consist of possibly incompressible
phases with finite material contrast. In section 2.1, we establish that the
velocity field in this type of layered material is phase-wise affine, if the
dissipation potentials of the phases are strictly rank-one convex. We use this
result and follow an analytical procedure described by Milton [69, §9] to
derive the closed form linear homogenization function of the considered
layered emulsions in section 2.2. To treat infinite material contrast, we
study a particular type of layered material, namely, coated layered materials
(CLMs), in section 2.3. We obtain closed form homogenization functions for
CLMs with kinematically unconstrained or incompressible phases. With the
objective of using CLMs as FDMN building blocks in the presence of rigid
fibers, we derive conditions for the required number of layering directions
and their orientation, such that the effective properties of CLM are non-
singular. We do so for CLMs with incompressible and rigid phases.

In section 3, we present the FDMN architecture as an extension of the (direct)
DMN architecture. An FDMN arises by replacing the rank-one laminates of a
(direct) DMN with rank-one layered materials capable of treating fluids, and
non-singular CLMs capable of treating infinite material contrast. Also, we
discuss the material sampling, the offline training procedure, as well as the
online evaluation of FDMNs in the context of incompressible rigid fiber
suspensions. We apply FDMNs to predict the viscous response of fiber polymer
suspensions in section 4, demonstrating both the ability to handle
incompressible fluid phases and infinite material contrast. We use FFT-based
computational homogenization techniques to generate linear training data and
non-linear validation data for shear-thinning fiber suspensions with a Cross-
type matrix behavior, and a fiber volume fraction of 25%. We consider 31
different fiber orientation states and a variety of incompressible
elongational and shear flows. Using the computational data, we train FDMNs for
all fiber orientation states, and study the prediction accuracy of the FDMNs
in the non-linear case. For all investigated load cases and microstructures,
the FDMNs achieve validation errors below 4.31% over a wide range of shear
rates relevant to engineering processes. Not considering the time required to
generate the training data and to train an FDMN, the FDMNs achieve speedup
factors between 11785 and 17225 compared to FFT-based computations. Finally,
we compare the prediction accuracy of the FDMN based approach with a different
machine learning approach by Sterr et al. [51], and find that FDMNs offer
better accuracy and flexibility at a higher computational cost for the
considered material and flow scenarios.

##  2 Homogenization of layered emulsions with infinite material contrast

###  2.1 Phase-wise affine displacement and velocity fields in laminates and
emulsions

Like for (direct) DMNs, the building blocks of an FDMN should admit linear
homogenization functions in closed form, and there should exist efficient
solution schemes to compute the stress response of the building blocks in case
their phases are non-linear. In analogy to the laminate based architecture of
the (direct) DMN, we thus consider layered emulsions as the building blocks of
the FDMN architecture. Like laminates, layered emulsions consist of multiple,
possibly non-linear, fluid materials, that are arranged such that every phase
boundary between the phases is orthogonal to the layering direction
$\mbox{\boldmath$n$}\in\mathcal{S}^{2}$, on the 2-sphere $\mathcal{S}^{2}$.
For the offline training and the online evaluation of an FDMN with layered
emulsions as building blocks, knowing the homogenization functions of layered
emulsions with linear and non-linear phases is necessary. Primarily, we look
for similarities between the homogenization of layered emulsions and solid
laminates that can be leveraged. In laminates, the displacement fields are
phase-wise affine, and we wish to know if the velocity fields in layered
emulsions belong to a particular class of fields as well. If the velocity
fields are phase-wise affine, established solution techniques on (direct) DMN
architectures [54] can be employed for the online evaluation of the FDMN, and
convenient linear homogenization equations can be employed for the offline
training. We follow a procedure detailed by Kabel et al. [70] on the
homogenization of laminates, and establish the existence and uniqueness of
phase-wise affine velocity fields in a particular type of emulsion. We
consider a $K$-phase layered emulsion with layering direction $n$, which,
similar to a laminate, is constructed by arranging $K$-phases so that the
direction $n$ is orthogonal to all phase boundaries. The emulsion occupies a
rectangular and periodic volume $Y\subset\mathds{R}^{3}$, and consists of non-
Newtonian phases with rank-one convex dissipation potential densities
${\Psi_{k}:(\mathds{R}^{3})^{\otimes 2}\rightarrow\mathds{R}}$, such that

$\Psi_{k}(\mbox{\boldmath$L$}+\beta\mbox{\boldmath$a$}\otimes\mbox{\boldmath$b$})\leq\beta\Psi_{k}(\mbox{\boldmath$L$}+\mbox{\boldmath$a$}\otimes\mbox{\boldmath$b$})+(1-\beta)\Psi_{k}(\mbox{\boldmath$L$})\quad\forall\,\beta\in[0,1],\,\mbox{\boldmath$a$},\mbox{\boldmath$b$}\in\mathds{R}^{3}\backslash\\{\mbox{\boldmath$0$}\\},$
(2.1)

where $\mbox{\boldmath$L$}\in(\mathds{R}^{3})^{\otimes 2}$ denotes the
velocity gradient, and the operator $(\cdot)^{\otimes a}$ constructs a tensor
space of $a$-th order. Also, $Y_{k}\subset Y$ denotes the sub volume of the
$k$-th phase, and the phases do not intersect, i.e.,

$Y_{j}\cap Y_{k}=\emptyset\quad\text{if}\quad j\neq k,$ (2.2)

and form the volume $Y$, such that

$Y=\bigcup_{j=1}^{K}Y_{j}.$ (2.3)

Let the operator $\left\langle{\cdot}\right\rangle_{Y}$ denote volume
averaging of a quantity over a volume $Y$, i.e.,

$\left\langle{\cdot}\right\rangle_{Y}\equiv\frac{1}{|Y|}\int_{Y}(\cdot)\mathop{}\\!\mathrm{d}Y(\mbox{\boldmath$x$})\quad\text{with}\quad\
|Y|\equiv\int_{Y}\mathop{}\\!\mathrm{d}Y(\mbox{\boldmath$x$}).$ (2.4)

Then, we express the velocity field
$\mbox{\boldmath$v$}:Y\rightarrow\mathds{R}^{3}$ as

$\mbox{\boldmath$v$}=\bar{\mbox{\boldmath$L$}}\mbox{\boldmath$x$}+\hat{\mbox{\boldmath$v$}},\quad\text{with}\quad\left\langle{\nabla\hat{\mbox{\boldmath$v$}}}\right\rangle_{Y}=\mbox{\boldmath$0$},$
(2.5)

where $\bar{\mbox{\boldmath$L$}}\in(\mathds{R}^{3})^{\otimes 2}$ stands for
the prescribed effective velocity gradient and $\hat{\mbox{\boldmath$v$}}$ is
the velocity fluctuation field. Accordingly, the periodic dissipation
potential density function ${\Psi:Y\times(\mathds{R}^{3})^{\otimes
2}\rightarrow\mathds{R}}$ of the emulsion is defined through

$\Psi(\mbox{\boldmath$x$},\nabla^{\sf{s}}\mbox{\boldmath$v$})=\sum_{i}^{K}\chi_{k}(\mbox{\boldmath$x$})\Psi_{k}(\mbox{\boldmath$x$},\nabla^{\sf{s}}\mbox{\boldmath$v$}),$
(2.6)

where $\chi_{k}:Y\rightarrow\\{0,1\\}$ refers to the characteristic function
of the sub volume $Y_{k}$, and $\nabla^{\sf{s}}$ stands for the symmetrized
gradient. Let the flow inside the emulsion be governed by Stokes flow, i.e.,
the advectorial forces in the layered emulsion are small in contrast to the
viscous forces, then the steady state balance of mass

${\mathrm{div}\ }{(\rho\mbox{\boldmath$v$})}=0,$ (2.7)

where $\rho:Y:\rightarrow\mathds{R}$ denotes the mass density field, is
satisfied. Additionally, the stress field

$\mbox{\boldmath$\sigma$}=\mathcal{S}(\cdot,\bar{\mbox{\boldmath$D$}}+\nabla^{s}\hat{\mbox{\boldmath$v$}})$
(2.8)

in terms of the effective strain rate tensor
$\bar{\mbox{\boldmath$D$}}=(\bar{\mbox{\boldmath$L$}}+{\bar{\mbox{\boldmath$L$}}}^{\sf{T}})/2$,
and the associated stress operator $\mathcal{S}$ on the space of symmetric
second order tensors $\sf{Sym}(3)\subset(\mathds{R}^{3})^{\otimes 2}$ [71]

$\displaystyle\mathcal{S}:Y\times\sf{Sym}(3)$
$\displaystyle\rightarrow\sf{Sym}(3),$ (2.9)
$\displaystyle(\mbox{\boldmath$x$},\mbox{\boldmath$D$})$
$\displaystyle\mapsto{\partialderivative{\Psi}{\mbox{\boldmath$D$}}}{(\mbox{\boldmath$x$},\mbox{\boldmath$D$})},$
(2.10)

satisfies the balance of linear momentum without inertial effects and volume
forces

${\mathrm{div}\ }{\mbox{\boldmath$\sigma$}}=\mbox{\boldmath$0$}.$ (2.11)

Therefore, the effective response of the emulsion is obtained from the
minimization problem

$\left\langle{\Psi(\mbox{\boldmath$x$},\bar{\mbox{\boldmath$D$}}+\nabla^{\sf{s}}\hat{\mbox{\boldmath$v$}})}\right\rangle_{Y}\longrightarrow\min_{\hat{\bm{v}}\in\mathcal{A}},$
(2.12)

where $\mathcal{A}$ stands for the admissible set

$\mathcal{A}=\left\\{\hat{\mbox{\boldmath$v$}}:Y\rightarrow\mathds{R}^{3}\quad\middle|\quad{\mathrm{div}\
}{(\rho\mbox{\boldmath$v$})}=0\quad\text{and}\quad\left\langle{\nabla\hat{\mbox{\boldmath$v$}}}\right\rangle_{Y}=\mbox{\boldmath$0$}\right\\}$
(2.13)

Hence, assuming there is no slip between the phases, the proof by Kabel et al.
[70, §2] extends to layered emulsions with finite material contrast.
Therefore, if the dissipation potential densities $\Psi_{k}$ are strictly
rank-one convex, there exists a unique minimizer $\hat{\mbox{\boldmath$v$}}$
in the class of phase-wise affine velocity fields. As a special case, this is
also true for incompressible Stokes flows, for which the balance of mass
simplifies to

${\mathrm{div}\ }{\mbox{\boldmath$v$}}=0.$ (2.14)

In the following, we use the term layered emulsion only for the type of
emulsion defined above.

###  2.2 Two-phase layered emulsions with Newtonian phases

Similar to two-phase laminates of linear elastic solid materials, we wish to
find a convenient closed form expression for the linear homogenization
function of two-phase layered emulsions. With the goal of using layered
emulsions as building blocks of an FDMN, we dedicate this section to deriving
such a closed form expression for two-phase layered emulsions with
kinematically unconstrained or incompressible Newtonian phases. For this
purpose, we employ the results of the previous section 2.1, and extend an
existing homogenization approach for two-phase laminates of linear elastic
phases by Milton [69, §9]. To emphasize the link between the linear
homogenization of laminates and layered emulsions we define the class of
linear materials $\mathcal{M}$, and study layered materials of solids and
fluids in parallel. A material $m_{k}$ in the class of linear materials
$\mathcal{M}$ is characterized by the primal material tensor
${\mathbb{M}}_{k}$ or its dual tensor
${\mathbb{K}}_{k}\equiv{\mathbb{M}}_{k}^{-1}$, such that

$({\mathbb{M}}_{k},{\mathbb{K}}_{k})\in\\{(\sf{M}^{+})^{2},(\sf{M}_{0}^{+})^{2}\\}.$
(2.15)

Here, the convex cones $\sf{M}^{+}\subset(\mathds{R}^{3})^{\otimes 4}$ and
$\sf{M}_{0}^{+}\subset(\mathds{R}^{3})^{\otimes 4}$ comprise all fourth order
tensors $\mbox{\boldmath$X$}\in(\mathds{R}^{3})^{\otimes 4}$, which are
positive definite on the vector spaces of symmetric second order tensors
$\sf{Sym}(3)$ and traceless symmetric second order tensors
${\sf{Sym}_{0}(3)}\subset\sf{Sym}(3)$, respectively. Additionally, a tensor
${\mathbb{X}}\in\\{\sf{M}^{+},\sf{M}_{0}^{+}\\}$ has minor and major
symmetries

${\mathbb{X}}={\mathbb{X}}^{\sf{T_{H}}}={\mathbb{X}}^{\sf{T_{L}}}={\mathbb{X}}^{\sf{T_{R}}},$
(2.16)

and a tensor ${\mathbb{X}}\in\sf{M}_{0}^{+}$ maps the second order unit tensor
$I$ to zero, such that

$\ker{\mathbb{X}}=\mbox{\boldmath$I$},\quad\text{i.e.,}\quad{\mathbb{X}}[\mbox{\boldmath$I$}]=\mbox{\boldmath$0$}.$
(2.17)

Here, we define the map ${\mathbb{T}}[\mbox{\boldmath$Y$}]$ in component form
and in the standard basis of $\mathds{R}^{3}$ as

${\mathbb{T}}[\mbox{\boldmath$Y$}]\mathrel{\stackon[1.5pt]{=}{\stretchto{\scalerel*[width("=")]{\wedge}{\rule{4.30554pt}{12.91663pt}}}{0.5ex}}}T_{ijkl}Y_{kl},\quad\text{for
some}\quad\mbox{\boldmath$Y$}\in(\mathds{R}^{3})^{\otimes
2},{\mathbb{T}}\in(\mathds{R}^{3})^{\otimes 4}.$ (2.18)

Any material $m_{k}\in\mathcal{M}$ follows the linear constitutive equations

$\mbox{\boldmath$\sigma$}^{co}={\mathbb{M}}_{k}[\nabla^{\sf{s}}\mbox{\boldmath$w$}]\quad\text{and}\quad\nabla^{\sf{s}}\mbox{\boldmath$w$}={\mathbb{K}}_{k}[\mbox{\boldmath$\sigma$}^{co}],$
(2.19)

if the material $m_{k}$ is not subjected to kinematic constraints, or the
constitutive equation

$\mbox{\boldmath$\sigma$}^{ico}=\mbox{\boldmath$\tau$}-p\mbox{\boldmath$I$}={\mathbb{M}}_{k}[\nabla^{\sf{s}}\mbox{\boldmath$w$}]-p\mbox{\boldmath$I$},\quad\text{and}\quad\nabla^{\sf{s}}\mbox{\boldmath$w$}={\mathbb{K}}_{k}[\mbox{\boldmath$\tau$}],$
(2.20)

if the material $m_{k}$ is incompressible. Here, $w$ denotes the field
associated with the type of material $m_{k}$, i.e., the displacement field $u$
for linear elastic solids or the velocity field $v$ for linearly viscous
fluids. For the constitutive equation (2.20) we use an additive split of the
stress field $\sigma$

$\mbox{\boldmath$\sigma$}=\mbox{\boldmath$\tau$}-p\mbox{\boldmath$I$},$ (2.21)

into the deviatoric stress field
$\mbox{\boldmath$\tau$}:Y\rightarrow{\sf{Sym}_{0}(3)}$ and the pressure field
$p:Y\rightarrow\mathds{R}$. Here, we use $Y\subset\mathds{R}^{3}$ to denote
the volume occupied by a two-phase layered material, i.e., a laminate or a
layered emulsion, and $Y_{1},Y_{2}\subset Y$ to denote the volume occupied by
the two phases $m_{1},m_{2}\in\mathcal{M}$ of the layered material. We define
the effective properties
$(\bar{{\mathbb{M}}},\bar{{\mathbb{K}}})\in\\{(\sf{M}^{+})^{2},(\sf{M}_{0}^{+})^{2}\\}$
of the layered material in the sense of the Hill-Mandel condition [72, 73]

$~{}\left\langle{\mbox{\boldmath$\sigma$}\cdot\nabla^{\sf{s}}\mbox{\boldmath$w$}}\right\rangle_{Y}=\left\langle{\mbox{\boldmath$\sigma$}}\right\rangle_{Y}\cdot\left\langle{\nabla^{\sf{s}}\mbox{\boldmath$w$}}\right\rangle_{Y}\equiv\bar{{\mathbb{M}}}\cdot\left(\left\langle{\nabla^{\sf{s}}\mbox{\boldmath$w$}}\right\rangle_{Y}\otimes\left\langle{\nabla^{\sf{s}}\mbox{\boldmath$w$}}\right\rangle_{Y}\right).$
(2.22)

Then, for perfectly bonded laminates of linear solid materials without
kinematic constraints, a closed form expression for the primal or dual
effective properties of the laminate
$\bar{{\mathbb{A}}}\in\\{\bar{{\mathbb{M}}},\bar{{\mathbb{K}}}\\}$ can be
derived using the procedure detailed by Milton [69, §9], such that

$\left(\left(\bar{{\mathbb{A}}}-{\mathbb{A}}_{0}\right)^{-1}+{\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{A}}_{0}}(\mbox{\boldmath$n$})\right)^{-1}=\left\langle\left(\left({\mathbb{A}}(\mbox{\boldmath$x$})-{\mathbb{A}}_{0}\right)^{-1}+{\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{A}}_{0}}(\mbox{\boldmath$n$})\right)^{-1}\right\rangle.$
(2.23)

Here, the inverse is taken on the space $\sf{M}^{+}$ in which
$\bar{{\mathbb{A}}}\in\sf{M}^{+}$ lives, ${\mathbb{A}}:Y\rightarrow\sf{M}^{+}$
denotes the field of the phase wise primal or dual material properties

${\mathbb{A}}:Y\rightarrow\sf{M}^{+},\quad\mbox{\boldmath$x$}\mapsto\begin{cases}{\mathbb{A}}_{1}&\text{if
$\mbox{\boldmath$x$}\in Y_{1}$},\\\ {\mathbb{A}}_{2}&\text{if
$\mbox{\boldmath$x$}\in Y_{2}$},\end{cases}$ (2.24)

and ${\mathbb{A}}_{0}\in\sf{M}^{+}$ refers to an arbitrary reference material
tensor such that the term
$({\mathbb{A}}(\mbox{\boldmath$x$})-{\mathbb{A}}_{0})$ is invertible. Also,
the the operator
${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{A}}_{0}}(\mbox{\boldmath$n$})\in(\mathds{R}^{3})^{\otimes
4}$ is defined as

${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{A}}_{0}}(\mbox{\boldmath$n$})={\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{I}}}(\mbox{\boldmath$n$})\left[{\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{I}}}(\mbox{\boldmath$n$}){\mathbb{A}}_{0}{\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{I}}}(\mbox{\boldmath$n$})\right]^{\dagger}{\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{I}}}(\mbox{\boldmath$n$}),$
(2.25)

where $(\cdot)^{\dagger}$ stands for the Moore–Penrose pseudoinverse, and the
operators
${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{A}}_{0}}(\mbox{\boldmath$n$})$
and
${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{I}}}(\mbox{\boldmath$n$})$
project onto the same subspace with respect to different inner products. More
specifically, the operator
${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{A}}_{0}}(\mbox{\boldmath$n$})$
is an orthogonal projector with respect to the ${\mathbb{A}}_{0}$-weighted
inner product

$\langle\mbox{\boldmath$X$},\mbox{\boldmath$X$}\rangle_{{\mathbb{A}}_{0}}=\mbox{\boldmath$X$}\cdot{\mathbb{A}}_{0}[\mbox{\boldmath$X$}],\quad\mbox{\boldmath$X$}\in\sf{Sym}(3),$
(2.26)

and the operator
${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{I}}}(\mbox{\boldmath$n$})$
is an orthogonal projector with respect to the Frobenius inner product

$\langle\mbox{\boldmath$X$},\mbox{\boldmath$X$}\rangle_{\sf{F}}=\mbox{\boldmath$X$}\cdot{\mathbb{I}}[\mbox{\boldmath$X$}]=\mbox{\boldmath$X$}\cdot\mbox{\boldmath$X$},\quad\mbox{\boldmath$X$}\in\sf{Sym}(3),$
(2.27)

where ${\mathbb{I}}\in(\mathds{R}^{3})^{\otimes 4}$ denotes the fourth order
identity tensor. The operator
${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{I}}}(\mbox{\boldmath$n$})$
encodes restrictions imposed by the periodicity of the laminate, the perfect
bonding of the phases, the momentum balance (2.11), and the constitutive
equation (2.19). Depending on whether we wish to compute the effective primal
or the effective dual properties of the laminate, a different operator
${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{I}}}(\mbox{\boldmath$n$})$
is required. For the homogenization of the primal effective properties, i.e.,
$\bar{{\mathbb{A}}}=\bar{{\mathbb{M}}}$, the operator
${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{I}}}(\mbox{\boldmath$n$})$
equals the operator
${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{1}^{\sf{co}}(\mbox{\boldmath$n$})$.
For the homogenization of the dual effective properties, i.e.,
$\bar{{\mathbb{A}}}=\bar{{\mathbb{K}}}$, the operator
${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{I}}}(\mbox{\boldmath$n$})$
equals the operator
${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{2}^{\sf{co}}(\mbox{\boldmath$n$})$.
The operators
${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{1}^{\sf{co}}$
and
${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{2}^{\sf{co}}$,
are defined through their action on a second order tensor
$\mbox{\boldmath$X$}\in\sf{Sym}(3)$

${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{1}^{\sf{co}}(\mbox{\boldmath$n$})[\mbox{\boldmath$X$}]=2(\mbox{\boldmath$X$}\mbox{\boldmath$n$})\otimes_{\sf{s}}\mbox{\boldmath$n$}-(\mbox{\boldmath$n$}\cdot\mbox{\boldmath$X$}\mbox{\boldmath$n$})\mbox{\boldmath$n$}\otimes\mbox{\boldmath$n$}\quad\text{and}\quad{\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{2}^{\sf{co}}(\mbox{\boldmath$n$})[\mbox{\boldmath$X$}]=\mbox{\boldmath$X$}-{\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{1}^{\sf{co}}(\mbox{\boldmath$n$})\mbox{\boldmath$X$},$
(2.28)

where $\otimes_{\sf{s}}$ denotes the symmetrized dyadic product. The operators
${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{1}^{\sf{co}}$
and
${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{2}^{\sf{co}}$
project onto the subspaces

$\mathcal{E}^{\sf{co}}=\left\\{\sf{Sym}(3)\ni\mbox{\boldmath$X$}=\mbox{\boldmath$d$}\otimes\mbox{\boldmath$n$}+\mbox{\boldmath$n$}\otimes\mbox{\boldmath$d$}\quad\middle|\quad\mbox{\boldmath$d$}\in\mathds{R}^{3}\right\\}\quad\text{and}\quad\mathcal{S}^{\sf{co}}=\left\\{\sf{Sym}(3)\ni\mbox{\boldmath$X$}\quad\middle|\quad\mbox{\boldmath$X$}\mbox{\boldmath$n$}=\mbox{\boldmath$0$}\right\\},$
(2.29)

and arise from the fact that the jumps of the field
${\llbracket{\nabla^{\sf{s}}\mbox{\boldmath$w$}}\rrbracket}\in\mathcal{E}^{\sf{co}}$

${\llbracket{\nabla^{\sf{s}}\mbox{\boldmath$w$}}\rrbracket}={\llbracket{\nabla\mbox{\boldmath$w$}\,\mbox{\boldmath$n$}}\rrbracket}\otimes_{\sf{s}}\mbox{\boldmath$n$},$
(2.30)

and the jump of the stress
${\llbracket{\mbox{\boldmath$\sigma$}}\rrbracket}\in\mathcal{S}^{\sf{co}}$
live on mutually orthogonal subspaces [69, §9], because it holds that

${\llbracket{\mbox{\boldmath$\sigma$}\mbox{\boldmath$n$}}\rrbracket}=\mbox{\boldmath$0$}.$
(2.31)

Here, we define the jump ${\llbracket{q}\rrbracket}$ of a quantity $q$ across
a phase boundary with normal vector $n$ as

${\llbracket{q}\rrbracket}=q^{+}-q^{-},$ (2.32)

where $q^{+}$ and $q^{-}$ denote the outer and inner limits of the quantity
$q$ with regard to the outside facing normal vector $n$, respectively.

Interestingly, the constitutive equations for kinematically unconstrained
linear elastic solids and linearly viscous fluids share the same structure,
see equation (2.19). Additionally, equations (2.30) and (2.31) are satisfied
for layered emulsions by definition, because there is no slip between the
phases, and the momentum balance (2.11) is satisfied. Thus, the homogenization
problem for the considered layered emulsions has the same structure as the
homogenization problem for laminates. As a result, the effective properties
$\bar{{\mathbb{A}}}\in\sf{M}^{+}$ of a two-phase layered emulsion of
kinematically unconstrained materials $m_{1},m_{2}\in\mathcal{M}$ can be
computed with equation (2.23). If the kinematic field $w$ is solenoidal such
that the incompressible balance of mass (2.14) is satisfied, it follows from
equation (2.30) that

${\llbracket{\nabla^{\sf{s}}\mbox{\boldmath$w$}\,\mbox{\boldmath$n$}}\rrbracket}\cdot\mbox{\boldmath$n$}=0.$
(2.33)

Also, if the materials $m_{1},m_{2}\in\mathcal{M}$ are incompressible, it
follows from equations (2.21) and (2.31) that

${\llbracket{\mbox{\boldmath$\sigma$}\mbox{\boldmath$n$}}\rrbracket}=\mbox{\boldmath$0$}\quad\iff\quad{\llbracket{\mbox{\boldmath$\tau$}\mbox{\boldmath$n$}}\rrbracket}={\llbracket{p}\rrbracket}\mbox{\boldmath$n$}.$
(2.34)

Therefore, the jumps
${\llbracket{\nabla^{\sf{s}}\mbox{\boldmath$w$}}\rrbracket}\in\mathcal{E}^{\sf{ico}}$
and
${\llbracket{\mbox{\boldmath$\sigma$}}\rrbracket}\in\mathcal{S}^{\sf{ico}}$
live in the mutually orthogonal spaces

$\mathcal{E}^{\sf{ico}}=\left\\{{\sf{Sym}_{0}(3)}\ni\mbox{\boldmath$X$}=2\,\mbox{\boldmath$d$}\otimes_{\sf{s}}\mbox{\boldmath$n$}-2(\mbox{\boldmath$d$}\cdot\mbox{\boldmath$n$})\mbox{\boldmath$n$}\otimes\mbox{\boldmath$n$}\quad\middle|\quad\mbox{\boldmath$d$}\in\mathds{R}^{3}\right\\}$
(2.35)

and

$\mathcal{S}^{\sf{ico}}=\left\\{{\sf{Sym}_{0}(3)}\ni\mbox{\boldmath$X$}\quad\middle|\quad\mbox{\boldmath$X$}\mbox{\boldmath$n$}=a\mbox{\boldmath$n$},\quad
a\in\mathds{R}\right\\}.$ (2.36)

The two operators
${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{1}^{\sf{ico}}(\mbox{\boldmath$n$})$
and
${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{2}^{\sf{ico}}(\mbox{\boldmath$n$})$,
defined through their action on a traceless second order tensor
${\mbox{\boldmath$X$}\in{\sf{Sym}_{0}(3)}}$,

${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{1}^{\sf{ico}}(\mbox{\boldmath$n$})[\mbox{\boldmath$X$}]=2(\mbox{\boldmath$X$}\mbox{\boldmath$n$})\otimes_{\sf{s}}\mbox{\boldmath$n$}-2(\mbox{\boldmath$n$}\cdot\mbox{\boldmath$X$}\mbox{\boldmath$n$})\mbox{\boldmath$n$}\otimes\mbox{\boldmath$n$}\quad\text{and}\quad{\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{2}^{\sf{ico}}(\mbox{\boldmath$n$})[\mbox{\boldmath$X$}]=\mbox{\boldmath$X$}-{\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{2}^{\sf{ico}}(\mbox{\boldmath$n$})[\mbox{\boldmath$X$}],$
(2.37)

project onto the spaces $\mathcal{E}^{\sf{ico}}$ and $\mathcal{S}^{\sf{ico}}$,
respectively. Consequently, the procedure detailed in Milton [69] can be used
to compute the effective properties $\bar{{\mathbb{A}}}\in\sf{M}_{0}^{+}$ of a
layered material with incompressible constituents. In this context the
homogenization equation (2.23) applies if the property field
${\mathbb{A}}:Y:\rightarrow\sf{M}_{0}^{+}$, the reference properties
${\mathbb{A}}_{0}\in\sf{M}_{0}^{+}$, and the operator
${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{I}}}(\mbox{\boldmath$n$})\in({\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{1}^{\sf{ico}}(\mbox{\boldmath$n$}),{\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{2}^{\sf{ico}}(\mbox{\boldmath$n$}))$
are selected based on whether the tensor $\bar{{\mathbb{A}}}$ represents
primal or dual effective properties. In summary, for a layered material
consisting of linear materials $m_{1},m_{2}\in\mathcal{M}$, the effective
material properties
$(\bar{{\mathbb{M}}},\bar{{\mathbb{K}}})\in\\{(\sf{M}^{+})^{2},(\sf{M}_{0}^{+})^{2}\\}$
can be computed with the homogenization equation (2.23). The choice of the
operator
${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{I}}}(\mbox{\boldmath$n$})$
depends on whether the phases are compressible or incompressible, and whether
the primal or dual effective tensors are sought.

###  2.3 Coated layered materials

Because we are interested in building FDMNs for suspensions of rigid
particles, the architecture of the FDMN should be capable of handling
infinitely viscous, i.e., rigid, phases. In this section, we leverage the
homogenization equation (2.23) to study the homogenization of repeatedly
layered materials in the presence of rigid phases. We define and use the term
rank of a layered material in the following to better distinguish between the
possible ways to construct layered materials. The rank $R\in\mathds{N}_{\geq
0}$ of a hierarchical, layered material defines the number of layering steps
that are needed to construct the material using layering directions
$\mbox{\boldmath$n$}_{r}\in\mathcal{S}^{2}$, such that the index ${r}$ starts
at one and ends at $R$, i.e., ${r}\in\\{1,...,R\\}$. By definition, a rank-$0$
layered material is simply a single phase material without any layering. Also,
the materials that are layered at each layering step need not be the same. For
example, let the first layer in a rank-2 layered material consist of the
materials $m_{1}$ and $m_{2}$. Then, the second layer may consist of the
materials $m_{1}$ and $m_{2}$ in the form of the first layer, and a third
material $m_{3}$.

Rank-one two-phase layered materials are unsuitable as DMN building blocks in
the case of infinite material contrast, because their effective properties are
always singular if one phase is singular. In search of a building block
capable of handling rigid phases, we study the effective properties of a
particular type of layered material following the terminology and the
analytical approach detailed by Milton [69, §9]. We refer to this particular
type of layered materials as coated layered materials (CLMs). To construct a
CLM, one material $m_{2}$, called the coating, is repeatedly layered onto
another material $m_{1}$, called the core. The layering process is divided
into multiple steps, such that the $r$-th layering step is applied in
direction $\mbox{\boldmath$n$}_{r}$, and the total number of layering steps
equals the rank $R$ of the CLM.

(a) Orthogonal layering directions

(b) Non-orthogonal layering directions

Figure 1: Rank-three, coated layered materials with orthogonal layering
directions (a), and non-orthogonal layering directions (b). The core material
is shown in red and the coating material is shown in gray. We hint at the
separation of length scales between the material layers through our choice of
the layer thicknesses.

As an example, two rank-three CLMs are depicted in Figure 1, where the
layering directions $\mbox{\boldmath$n$}_{r}\in\mathcal{N}_{1}$ of the
material shown in Figure 1(a) are

$\mathcal{N}_{1}=\left\\{(1,0,0)_{\\{\mbox{\boldmath$e$}_{i}\\}},(0,0,1)_{\\{\mbox{\boldmath$e$}_{i}\\}},(0,1,0)_{\\{\mbox{\boldmath$e$}_{i}\\}}\right\\},$
(2.38)

and the layering directions $\mbox{\boldmath$n$}_{r}\in\mathcal{N}_{2}$ of the
material shown in Figure 1(b) are

$\mathcal{N}_{2}=\left\\{(1,0,0)_{\\{\mbox{\boldmath$e$}_{i}\\}},\frac{1}{\sqrt{2}}\left(1,1,0\right)_{\\{\mbox{\boldmath$e$}_{i}\\}},\frac{1}{\sqrt{2}}\left(1,0,1\right)_{\\{\mbox{\boldmath$e$}_{i}\\}}\right\\}.$
(2.39)

As described, the gray coating material $m_{2}$ is repeatedly applied in
layers, such that for ${r}>1$, the $r$-th sub layer of the rank-three CLM is
constructed by layering a rank-(${r}-1$) layered material with the coating
material $m_{2}$. By definition, the $0$-th sub layer of the rank-three CLM is
purely made of the core material $m_{1}$. We wish to compute the effective
properties $\bar{{\mathbb{A}}}$ of a coated layered material by repeatedly
applying the homogenization equation (2.23). For the recursive application of
the rank-one homogenization formula (2.23) to be valid, the phases at each
homogenization step must be homogeneous on the length scale of the respective
layering step. This requires a separation of length scales, in the sense that
the thickness $\delta_{r}$ of the material sections in the rank-${r}$ sub
layer is drastically smaller than the thickness $\delta_{{r}+1}$ of the
material sections in the next sub layer, i.e., $\delta_{r}\ll\delta_{r+1}$. In
Figure 1, we hint at this separation of length scales between the material
layers through our choice of the layer thicknesses. Nonetheless, we retain a
visible degree of inhomogeneity, so that the layer structures can be
distinguished for illustration purposes.

Before we repeatedly apply equation (2.23) to derive an expression for the
effective properties $\bar{{\mathbb{A}}}$, we first simplify the expression
(2.23) for the case at hand. We let the reference material properties
${\mathbb{A}}_{0}$ in equation (2.23) approach the properties
${\mathbb{A}}_{2}$ of the coating material $m_{2}$, resulting in the equation
[69, §9]

$(1-f_{2})(\bar{{\mathbb{A}}}-{\mathbb{A}}_{2})^{-1}=\left({\mathbb{A}}_{1}-{\mathbb{A}}_{2}\right)^{-1}+f_{2}{\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{A}}_{2}}(\mbox{\boldmath$n$}).$
(2.40)

Here, the coating material $m_{2}$ must not be rigid or a void, the tensor
$({\mathbb{A}}_{1}-{\mathbb{A}}_{2})$ must be invertible, and $f_{2}$ denotes
the volume fraction of the coating material $m_{2}$. Also, the fourth order
tensor
${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{A}}_{2}}(\mbox{\boldmath$n$}_{r})\in(\mathds{R}^{3})^{\otimes
4}$ is defined by equation (2.25) as

${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{A}}_{2}}(\mbox{\boldmath$n$}_{r})={\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{I}}}(\mbox{\boldmath$n$}_{r})\left[{\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{I}}}(\mbox{\boldmath$n$}_{r}){\mathbb{A}}_{2}{\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{I}}}(\mbox{\boldmath$n$}_{r})\right]^{\dagger}{\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{I}}}(\mbox{\boldmath$n$}_{r}),$
(2.41)

where the inverse is taken on the subspace onto which
${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{I}}}$
projects. To compute the effective properties $\bar{{\mathbb{A}}}$ of a
rank-$R$ CLM, we recursively apply equation (2.40) for the effective material
properties of rank-one layered materials, such that [69, §9]

$(1-f_{2})(\bar{{\mathbb{A}}}-{\mathbb{A}}_{2})^{-1}=({\mathbb{A}}_{1}-{\mathbb{A}}_{2})^{-1}+f_{2}\sum_{r=1}^{R}c_{r}{\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{A}}_{2}}(\mbox{\boldmath$n$}_{r}).$
(2.42)

Here, the scalars $c_{r}$ are computed from the volume fractions $f^{(r)}_{1}$
of the core material $m_{1}$ in the rank-${r}$ sub layer of the rank-three
CLM, such that

$c_{r}=\frac{f^{({r}-1)}_{1}-f^{({r})}_{1}}{f_{2}},$ (2.43)

i.e., the properties

$\sum_{r=1}^{R}c_{r}=1\quad\text{and}\quad c_{r}>0,$ (2.44)

hold. Because the rank-$0$ sub layer consists solely of the core material
$m_{1}$ by definition, the volume fraction $f^{(0)}=1$ is one. Equation (2.42)
is essential to compute the effective properties of layered materials with
infinite material contrast, and thus for the treatment of such materials using
a DMN architecture. In Appendix A, we employ equation (2.42) to study the
conditions under which the effective properties $\bar{{\mathbb{A}}}$ of a CLM
are singular when the core material $m_{1}$ is rigid. With the objective of
applying FDMNs to suspensions of rigid fibers in mind, we summarize the
results pertaining to rank-$R$ CLMs with a rigid core material $m_{1}$ and an
incompressible coating material $m_{2}$ with the following statements.

  1. 1.

The effective properties $\bar{{\mathbb{A}}}$ of a CLM of rank $R<3$ are
always singular.

  2. 2.

For rank-3 CLMs, the effective properties $\bar{{\mathbb{A}}}$ are singular if
at least one layering direction is orthogonal to two other layering
directions, or if at least two layering directions are collinear.

  3. 3.

Rank-$3$ CLMs with mutually non-orthogonal and mutually non-collinear layering
directions $\mbox{\boldmath$n$}_{1},\mbox{\boldmath$n$}_{2},$ and
$\mbox{\boldmath$n$}_{3}$ have non-singular effective properties. Thus, in the
context of this article, we require

$0<|\mbox{\boldmath$n$}_{1}\cdot\mbox{\boldmath$n$}_{2}|<1,\quad
0<|\mbox{\boldmath$n$}_{1}\cdot\mbox{\boldmath$n$}_{3}|<1,\quad
0<|\mbox{\boldmath$n$}_{2}\cdot\mbox{\boldmath$n$}_{3}|<1.$ (2.45)

The statements 1-3 are of critical importance for the treatment of
incompressible suspensions of rigid particles using a DMN architecture. It
follows that CLMs with an incompressible coating material and rigid core
material should be at least of rank $R=3$ and have non-collinear and mutually
non-orthogonal layering directions to be non-singular. Therefore, we leverage
such non-singular-singular rank-3 CLMs in the FDMN architecture to treat
suspensions of rigid fibers, as we explain in section 3.2.

##  3 Identifying deep material networks for suspensions of rigid fibers

###  3.1 Architecture of direct deep material networks

Because the FDMN architecture extends the (direct) DMN framework, we briefly
summarize the direct DMN architecture in this section, before presenting the
FDMN architecture in the next section 3.2. For in depth discussions and
explanations of the direct DMN architecture, we refer the reader to articles
by Gajek et al. [54, 46, 59]. A two-phase direct DMN is defined as a perfect
binary tree, consisting of two-phase laminates $\mathcal{B}^{i}_{k}$ as nodes,
see Figure 2(a). The laminates $\mathcal{B}^{i}_{k}$ are defined by a single
lamination direction $\mbox{\boldmath$n$}_{k}^{i}$ and two volume fractions
$c_{k,1}^{1}$ and $c_{k,1}^{2}$. The DMN tree is ordered, rooted and has depth
$K$. We use the letter ${k}=1,...,K$ to denote the depth of a node and the
letter $i=1,...,2^{k-1}$ to indicate the horizontal position of a node in the
respective layer.

(a) Weight propagation

(b) Stiffness propagation

Figure 2: Weight and stiffness propagation (from the bottom to the top) in a
two-phase direct DMN [46] of depth K = 3, by Gajek et al. [59, Fig. 1],
licensed under CC BY 4.0

We parametrize a direct DMN by the two vectors
$\underline{n}\in(\mathds{R}^{3})^{2^{K}-1}$

$\underline{n}=(\mbox{\boldmath$n$}_{K}^{1},....,\mbox{\boldmath$n$}_{K}^{2^{K-1}},\mbox{\boldmath$n$}_{K-1}^{1},...,\mbox{\boldmath$n$}_{K-1}^{2^{K-2}},...,\mbox{\boldmath$n$}_{1}^{1}),$
(3.1)

and $\underline{w}\in\mathds{R}^{2^{K}}$

$\underline{w}=(w_{K+1}^{1},...,w_{K+1}^{2^{K}}).$ (3.2)

Instead of the laminate volume fractions, the direct DMN is parametrized in
terms of the weights $w_{K+1}^{i}$ to improve numerical stability during the
training process [52, 53]. The $i$-th weight in the $k$-th layer is computed
from the pairwise summation of the weights in the previous layer $k+1$, see
Figure 2(a), such that

$w_{k}^{i}=w_{k+1}^{2i-1}+w_{k+1}^{2i},$ (3.3)

and the volume fractions $c_{k,1}^{i}$ and $c_{k,2}^{i}$ follow from

$c_{k,1}^{i}=\frac{w_{k+1}^{2i-1}}{w_{k+1}^{2i-1}+w_{k+1}^{2i}},\quad\text{and}\quad
c_{k,2}^{i}=1-c_{k,1}^{i}.$ (3.4)

Additionally, the weights $w_{K+1}^{i}$ need to be positive and sum to unity,
such that

$w_{K+1}^{i}\geq 0\quad\text{and}\quad\sum_{i}^{2^{K}}w_{K+1}^{i}=1.$ (3.5)

The vectors $\underline{n}$ and $\underline{w}$ fully parametrize a direct DMN
of depth $K$ and are identified by fitting the linear homogenization function
of the DMN to the linear homogenization function of the microstructure of
interest. The fitting process is carried out by using supervised machine
learning, and we refer to this process as offline training. To evaluate the
linear homogenization function of the DMN during the offline training, input
stiffness pairs are propagated through the network of homogenization blocks
$\mathcal{B}^{i}_{k}$, see Figure 2(b). Consequently, if the stiffness
${\mathbb{C}}_{k}^{i}$ of a homogenization block $\mathcal{B}^{i}_{k}$ is
singular, this singular stiffness can propagate through the whole DMN, which
is undesirable for the offline training. Once the parameter vectors
$\underline{n}$ and $\underline{w}$ are known, they remain unchanged and are
used in the online evaluation to predict the non-linear response of the
studied material. In sections 3.4, 4.2, and 3.5, we discuss the offline
training and the online evaluation for the proposed FDMN architecture.

###  3.2 Architecture of flexible deep material networks

With the results of the previous sections 2.1-2.3 at hand, we present an
extended DMN architecture that is able to treat fluid phases and infinite
material contrast. Compared to the architecture of (direct) DMNs shown in
Figure 2, the lowest layer of laminates is replaced by a layer of rank-$R$
CLMs $\mathcal{C}_{R}$, see Figure 3.

(a) Weight propagation

(b) Material tensor propagation

Figure 3: Weight and material property propagation (from the bottom to the
top) in a two-phase FDMN [46] of depth $K=3$ with rank-$3$ CLMs

The rank $R$ of the CLM $\mathcal{C}_{R}$ and the restrictions on its layering
directions $\mbox{\boldmath$n$}_{r}$ are problem dependent, and chosen such
that the non-singularity condition (5.11) is satisfied. This guarantees that
only non-singular material tensors ${\mathbb{A}}_{k}^{i}$ are propagated
through the network. If all phases are kinematically unconstrained, solid, and
non-singular, rank-$1$ CLMs can be employed to recover the direct DMN
architecture. However, in the case of incompressible suspensions of rigid
particles, the CLMs are of rank $R=3$ and have mutually non-collinear and non-
orthogonal layering directions $\mbox{\boldmath$n$}_{r}$, as per statements
1-3. Additionally, the homogenization blocks $\mathcal{H}^{i}_{k}$ implement
the general homogenization function defined by equation (2.23) in which the
projector function
${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{I}}}$
depends on the physics of the investigated problem. This allows the treatment
of composites consisting of linear materials $m\in\mathcal{M}$, including
incompressible linearly viscous fluids. For problems involving two-phase
materials, $R+1$ weights and two materials are assigned to each CLM. The
number of layering directions $N_{\sf{n}}$ and the number of input weights
$N_{\sf{w}}$ in a two-phase FDMN of depth $K$ with rank-$R$ CLMs are

$~{}N_{\sf{n}}=2^{K-1}(R+1)-1,\quad N_{\sf{w}}=2^{K-1}(R+1).$ (3.6)

As an example, Figures 3(a) and 3(b) show the weight and material property
propagation of a two-phase FDMN of depth $K=3$ with two-phase CLMs of rank
$R=3$ as bottom layer. Of the four weights of a rank-$3$ CLM, one is
associated with the core material $m_{1}$, and three are associated with the
coating material $m_{2}$. This doubles the number of weights for an FDMN with
rank-$3$ CLMs as compared to a DMN of the same depth $K$. Additionally, each
rank-$3$ CLM block $\mathcal{C}_{3}$ has three layering directions
$\mbox{\boldmath$n$}_{r}$ instead of one layering direction like rank-one
laminates. To reduce the number of free parameters and ensure the CLMs are
non-singular, the relative angles between the layering directions
$\mbox{\boldmath$n$}_{r}$ are fixed, while their joint orientation in space is
determined through the offline training, as we explain in section 3.4 in more
detail. Because CLMs are based on the repeated layering of rank-one layered
materials, an FDMN of depth $K$ can be considered to be a special case of a
(direct) DMN of depth $K+(R-1)$ with generalized homogenization blocks. Then,
a CLM represents a subtree of the DMN of depth $K+(R-1)$, in which the weights
$w_{k}^{i}$ of certain phases are set to zero, such that a CLM emerges, see
Figure 4.

(a) Weight propagation

(b) Material tensor propagation

Figure 4: Weight and material property propagation (from the bottom to the
top) in a two-phase DMN of depth $K=3$, representing a two-phase rank-$3$ CLM

In particular, except for one phase, the weights associated with the singular
core material $m_{1}$ are set to zero. Additionally, $R$ weights associated
with the non-singular coating material $m_{2}$ must be non-zero, such that the
rank $R$ of the CLM is maintained. Technically, (direct) DMNs are therefore
able to treat infinite material contrast intrinsically. However, the weights
$w_{k}^{i}$ need to be identified during the offline training, and might not
always have appropriate values to generate CLMs with the required rank, such
that singular material tensors ${\mathbb{A}}_{k}^{i}$ might propagate through
the whole DMN. Also, if the layering normals are not chosen appropriately, the
singular material properties may propagate through the whole DMN as well.
Controlling the weights and normals via the implementation of the CLM blocks
instead of learning them is therefore preferable to prevent singularities
during the offline training. Furthermore, FDMNs of depth $K$ have fewer free
parameters compared to equivalent DMNs of depth $K+R-1$, which aids in the
offline training as well.

Hypothetically, the presented FDMN architecture could be applied to problems
with more than two phases as well. Let $N_{R}$ denote the number of non-
singular phases and $N_{S}$ denote the number of singular phases. Then, for
each singular phase, the FDMN needs to contain at least one CLM to which the
phase is assigned, and which has non-singular effective properties. This
requires a minimum of $N_{S}$ CLMs. Additionally, the coatings of the CLMs can
be constructed by mixing the non-singular phases, possibly by constructing
$N_{R}$-phase layered materials. However, the application of FDMNs to
materials with voids or more than two phases is beyond the scope of this
article, because we focus on suspensions of rigid fibers.

###  3.3 Material sampling

For materials with solid phases, DMNs are trained by sampling the linear
elastic homogenization function of the microstructure of interest with
randomly generated input materials [52, 53, 46, 59], and fitting the DMN
parameters to replicate the desired homogenization function. Similarly, in the
case of rigid particles suspended in viscous media, we train FDMNs to learn
the linear viscous homogenization function with the goal of predicting non-
linear viscous behavior. Consequently, we sample the linear viscous
homogenization function by using randomly generated, linearly viscous input
materials. We base our sampling strategy on the approach used by Gajek et al.
[46]. In the article by Gajek et al. [46] on glass fiber reinforced polyamide
composites, the authors chose the stiffness samples for the fiber material to
be isotropic, and the stiffness samples for to the polyamide matrix to be
transversely isotropic. This strategy leads to a smaller sample space compared
to sampling orthotropic stiffnesses for both the fiber and the matrix
materials, as suggested by Liu et al. [52, 53]. To reduce the computational
cost of sampling even further for the considered case of fiber suspension, we
only generate sample viscosities for the matrix material and consider the
suspended fibers to be rigid in all computations. To generate samples for the
matrix material, we follow Gajek et al. [46] and choose the sample viscosity
${\mathbb{V}}^{s}$ pertaining to the $s$-th sample to be transversely
isotropic, such that

${\mathbb{V}}^{s}=2\eta_{2}\left({\mathbb{P}}_{2}-g\mbox{\boldmath$A$}^{\prime}(\underline{a})\otimes\mbox{\boldmath$A$}^{\prime}(\underline{a})\right),$
(3.7)

where the scalar $g\in[0,1)$ defines the magnitude of the perturbation by the
unitary and deviatoric tensor
$\mbox{\boldmath$A$}^{\prime}(\underline{a})\in{\sf{Sym}_{0}(3)}$, the vector
$\underline{a}\in\mathds{R}^{4}$ collects four angles, and we consider the
scalar viscosities

$\eta_{2}=10^{p}\mathrm{GPa},\quad p\in[-3,3].$ (3.8)

For a fiber reinforced polyamide, Gajek et al. [46] employed an ansatz of the
form (3.7) to sample the material tangents of a J2-elastoplastic matrix
material. Similarly, with the objective of applying FDMNs to fiber suspensions
with a shear-thinning polyamide matrix as detailed in section 4.1, we use
equation (3.7) to sample the material tangents of a Cross-type (4.4) matrix
material. The tensor function $\mbox{\boldmath$A$}^{\prime}$ is defined by an
eigenvalue decomposition [46]

$\displaystyle\mbox{\boldmath$A$}^{\prime}:\mathds{R}^{4}$
$\displaystyle\rightarrow{\sf{Sym}_{0}(3)},$ (3.9)
$\displaystyle{(a_{1},a_{2},a_{3},a_{4})}^{\sf{T}}$
$\displaystyle\mapsto\mbox{\boldmath$Q$}(a_{2},a_{3},a_{4})\mbox{\boldmath$A$}(a_{1}){\mbox{\boldmath$Q$}(a_{2},a_{3},a_{4})}^{\sf{T}},$
(3.10)

into the orthogonal tensor $\mbox{\boldmath$Q$}(a_{2},a_{3},a_{4})$ and a
diagonal tensor $\mbox{\boldmath$A$}(a_{1})$ in terms of the vector
${\underline{a}={(a_{1},a_{2},a_{3},a_{4})}^{\sf{T}}}$. The tensor
$\mbox{\boldmath$A$}(a_{1})$ is constructed from the angle $a_{1}$ via [46]

$A$
$\displaystyle=\frac{1}{\sqrt{{\frac{\xi_{1}+1}{2\,\xi_{1}+1}}}}\mbox{\boldmath$\mathrm{diag}$}\left(\xi_{2}\cos\left(a_{1}\right),\xi_{2}\sin{\left(a_{1}\right)},\frac{1}{\sqrt{2}}\right),$
(3.11) $\displaystyle\xi_{1}$
$\displaystyle=\cos\left(a_{1}\right)\sin\left(a_{1}\right),\quad\xi_{2}=\frac{-\sqrt{2}}{2\cos{\left(a_{1}\right)}+2\sin{\left(a_{1}\right)}},$
(3.12)

and the tensor $\mbox{\boldmath$Q$}(a_{2},a_{3},a_{4})$ represents a rotation
around the direction $n$, such that

$\displaystyle\mbox{\boldmath$Q$}(a_{2},a_{3},a_{4})$
$\displaystyle:{\mathbb{R}}^{3}\rightarrow{\mathbb{R}}^{3},\
\mbox{\boldmath$x$}\mapsto\cos\left(a_{2}\right)\mbox{\boldmath$x$}+\sin\left(a_{2}\right)\mbox{\boldmath$n$}\times\mbox{\boldmath$x$}+(1-\cos\left(a_{2}\right))(\mbox{\boldmath$n$}\cdot\mbox{\boldmath$x$})\mbox{\boldmath$n$},$
(3.13) $n$
$\displaystyle=\sin\left(a_{3}\right)\cos\left(a_{4}\right)\mbox{\boldmath$e$}_{1}+\sin\left(a_{3}\right)\sin\left(a_{4}\right)\mbox{\boldmath$e$}_{2}+\cos\left(a_{3}\right)\mbox{\boldmath$e$}_{3}.$
(3.14)

Here, the operator $\mathrm{diag}$ constructs a second order tensor in the
standard basis $\mbox{\boldmath$e$}_{i}$ of $\mathds{R}^{3}$, and the
components of the vector $\underline{a}$ live on the intervals

$a_{1}\in[0,2\pi],\quad a_{2}-\sin(a_{2})\in[0,\pi],\quad
a_{3}\in[0,\pi],\quad a_{4}\in[0,2\pi].$ (3.15)

We sample the space of the variables $p$ and $\underline{a}$ using the
scrambled Sobol sequence [74, 75], and employ an FFT-based computational
procedure to compute the resulting effective viscosities
$\bar{{\mathbb{V}}}^{s}$ [38, 17], see section (4).

###  3.4 Offline training

Using the material sampling procedure described in the previous section 3.3,
we wish to optimize the parameters of the FDMN with regard to the training
data set. For the case of rigid fibers suspended in viscous media, we consider
FDMNs with rank-$3$ CLMs as bottom layer, which satisfy the requirements
derived in section A.2. We collect all layering directions of the FDMN in a
vector

$\underline{n}=(\mbox{\boldmath$n$}_{K}^{1},...,\mbox{\boldmath$n$}_{K}^{3(2^{K-1})},\mbox{\boldmath$n$}_{K-1}^{1},...,\mbox{\boldmath$n$}_{K-1}^{2^{K-2}},...,\mbox{\boldmath$n$}_{2}^{1},\mbox{\boldmath$n$}_{2}^{2},\mbox{\boldmath$n$}_{1}^{1}),$
(3.16)

where we used the fact that the $k$-th layer contains $2^{k-1}$ normals if
$k<K$, and the $K$-th layer contains $3(2^{K-1})$ normals. Because the
homogenization blocks in the $K$-th layer are rank-three CLMs, the number of
normals is different from the other layers. Each normal
$\mbox{\boldmath$n$}_{k}^{i}$ in the $k$-th layer for $k<K$ is constructed
using two angles $a_{k}^{i}$ and $b_{k}^{i}$ via

$\mbox{\boldmath$n$}_{k}^{i}=\sin(b_{k}^{i})\cos(a_{k}^{i})\mbox{\boldmath$e$}_{1}+\sin(b_{k}^{i})\sin(a_{k}^{i})\mbox{\boldmath$e$}_{2}+\cos(b_{k}^{i})\mbox{\boldmath$e$}_{3}.$
(3.17)

However, for the normals $\mbox{\boldmath$n$}_{K}^{i}$ in the $K$-th layer,
the procedure is different. To avoid singular effective properties during the
evaluation of the FDMN in accordance with the discussion in sections 2.3 and
A.2, we ensure that the three layering normals of a rank-three CLM block are
not mutually non-orthogonal and mutually non-collinear as follows. First, for
each rank-three CLM block $\mathcal{C}_{3}^{i}$, we compute the components
$Q^{i}_{op}$ of the orthogonal tensor $\mbox{\boldmath$Q$}^{i}$ in the
standard basis $\\{\mbox{\boldmath$e$}_{i}\\}$ using three Euler angles
$c_{K}^{i},d_{K}^{i},e_{K}^{i}\in[0,2\pi]$ via

$Q^{i}_{op}=\left[\begin{array}[]{ccc}1&0&0\\\
0&\cos(c_{K}^{i})&-\sin(c_{K}^{i})\\\ 0&\sin(c_{K}^{i})&\cos(c_{K}^{i})\\\
\end{array}\right]\left[\begin{array}[]{ccc}\cos(d_{K}^{i})&0&\sin(d_{K}^{i})\\\
0&1&0\\\ -\sin(d_{K}^{i})&0&\cos(d_{K}^{i})\\\
\end{array}\right]\left[\begin{array}[]{ccc}\cos(e_{K}^{i})&-\sin(e_{K}^{i})&0\\\
\sin(e_{K}^{i})&\cos(e_{K}^{i})&0\\\ 0&0&1\\\ \end{array}\right].$ (3.18)

Then, we compute the layering directions
$\mbox{\boldmath$n$}_{K}^{3(i-1)+1},\mbox{\boldmath$n$}_{K}^{3(i-1)+2},\mbox{\boldmath$n$}_{K}^{3(i-1)+3}$
of the rank-three CLM $\mathcal{C}_{3}^{i}$ via

$\mbox{\boldmath$n$}_{K}^{3(i-1)+1}=\mbox{\boldmath$q$}^{i}_{1},\quad\mbox{\boldmath$n$}_{K}^{3(i-1)+2}=\frac{1}{\sqrt{2}}\left(\mbox{\boldmath$q$}^{i}_{1}+\mbox{\boldmath$q$}^{i}_{2}\right),\quad\mbox{\boldmath$n$}_{K}^{3(i-1)+3}=\frac{1}{\sqrt{2}}\left(\mbox{\boldmath$q$}^{i}_{1}+\mbox{\boldmath$q$}^{i}_{3}\right),$
(3.19)

where $\mbox{\boldmath$q$}^{i}_{1}$, $\mbox{\boldmath$q$}^{i}_{2}$, and
$\mbox{\boldmath$q$}^{i}_{3}$ are the first, second, and third column vectors
of the tensor $\mbox{\boldmath$Q$}^{i}$ in the standard basis
$\mbox{\boldmath$e$}_{i}$, respectively. This leads to a total of $N_{\alpha}$
angles with

$N_{\alpha}=5\left(2^{K-1}\right)-2.$ (3.20)

The angles between the layering directions are fixed for each CLM
$\mathcal{C}_{3}$, but the three layering directions are rotated together
during training. Other procedures than the one presented here are viable as
well, as long as the constructed vectors $\mbox{\boldmath$n$}_{K}^{3(i-1)+1}$,
$\mbox{\boldmath$n$}_{K}^{3(i-1)+2}$, and $\mbox{\boldmath$n$}_{K}^{3(i-1)+3}$
are mutually non-collinear and non-orthogonal. Alternatively, one could
parametrize each of the three layering directions
$\mbox{\boldmath$n$}_{K}^{3(i-1)+1}$, $\mbox{\boldmath$n$}_{K}^{3(i-1)+2}$,
and $\mbox{\boldmath$n$}_{K}^{3(i-1)+3}$ with two angles each, and identify
the angles independently. Then, the effective properties of the CLM
$\mathcal{C}_{3}^{i}$ are non-singular almost surely, because the probability
of any two of the three layering directions being collinear or one vector
being orthogonal to the other two is zero. However, parametrizing each of the
three layering directions with two angles requires twice the amount of angles
in the CLM layer of the FDMN as compared to the presented procedure.

To define the linear homogenization function of the FDMN, we collect all
$N_{\alpha}$ angles of the layering directions in the vector
$\underline{\alpha}$

$\underline{\alpha}=(c_{K}^{1},d_{K}^{1},e_{K}^{1},...,c_{K}^{3(2^{K-1})},d_{K}^{3(2^{K-1})},e_{K}^{3(2^{K-1})},a_{K-1}^{1},b_{K-1}^{1},...,a_{K-1}^{2^{K-2}},b_{K-1}^{2^{K-2}},...,a_{1}^{1},b_{1}^{1}),$
(3.21)

and collect all $N_{\sf{w}}$ input weights of the FDMN in the vector

$\underline{w}=(w_{K+1}^{1},...,w_{K+1}^{N_{\sf{w}}}).$ (3.22)

We follow Gajek et al. [54] and define the loss function
$\mathcal{L}:\mathds{R}^{\otimes N_{\alpha}}\times\mathds{R}^{\otimes
N_{\sf{w}}}\rightarrow\mathds{R}$ in terms of the linear homogenization
function
$\mathcal{DMN}_{\sf{L}}:\mathcal{A}\times\mathcal{A}\times{\mathds{R}}^{\otimes
N_{\alpha}}\times{\mathds{R}}^{\otimes N_{\sf{w}}}\rightarrow\mathcal{A}$ of
the FDMN as

$\mathcal{L}\left(\underline{\alpha},\underline{w}\right)=\frac{1}{N_{\sf{b}}}\sqrt[q]{\sum_{s=1}^{N_{\sf{b}}}\left(\frac{\norm{\bar{{\mathbb{A}}}^{s}-\mathcal{DMN}_{\sf{L}}\left({\mathbb{A}}^{s}_{1},{\mathbb{A}}^{s}_{2},\underline{\alpha},\underline{w}\right)}_{p}}{\norm{\bar{{\mathbb{A}}}^{s}}_{p}}\right)^{q}}+\lambda\left(\sum_{i=1}^{N_{\sf{w}}}w_{K+1}^{i}-1\right)^{2},$
(3.23)

where $N_{\sf{b}}$ is the number of samples,
${\mathbb{A}}_{1}^{s},{\mathbb{A}}_{2}^{s}\in\mathcal{A}$ denote the material
tensors of the phases of a sample, and $\bar{{\mathbb{A}}}^{s}\in\mathcal{A}$
denotes the effective material tensor of a sample. Additionally, the class of
material tensors $\mathcal{A}\in\\{\sf{M},\sf{M}_{0}\\}$ depends on the
considered problem, we follow Gajek et al. [54] and choose the constants
$p,q\in\mathds{N}_{\geq 0}$ as $p=1$ and $q=10$, and the operator
$\norm*{\cdot}_{1}$ refers to the $\ell^{1}$-norm of the components in Mandel
notation. The second summand in equation (3.23) encodes the simplex constraint
on the weights of the DMN

$\sum_{i=1}^{N_{\sf{w}}}w_{K+1}^{i}=1,\quad\text{where}\quad w_{K+1}^{i}\geq
0,$ (3.24)

in the form of a quadratic penalty term with the penalty parameter
$\lambda>0$. We implement the offline training in the Python programming
language and use the machine learning framework PyTorch [76] to identify the
vectors $\underline{\alpha}$ and $\underline{w}$. To do so, we leverage the
automatic differentiation capabilities of PyTorch to solve the minimization
problem

$\mathcal{L}\left(\underline{\alpha},\underline{w}\right)\rightarrow\min_{\underline{\alpha},\underline{w}},$
(3.25)

using the RAdam[77] optimizer. We use mini batches of $N_{\sf{b}}$ samples to
compute the parameter updates

$\underline{\alpha}_{q+1}=\underline{\alpha}_{q}-\kappa_{q,\alpha}{\partialderivative{\mathcal{L}}{\underline{\alpha}}}(\underline{\alpha},\underline{w}),\quad\underline{w}_{q+1}=\underline{w}_{q}-\kappa_{q,w}{\partialderivative{\mathcal{L}}{\underline{w}}}(\underline{\alpha},\underline{w})$
(3.26)

for every epoch $q$. Also, we follow previous work [52, 53, 54] and use
perfect binary trees without tree compression [52, 53, 46]

###  3.5 Online evaluation

In section 2.1, we showed that the velocity fields in layered emulsions are
phase-wise affine, and thus the online evaluation procedure of direct DMNs can
be employed for FDMNs as well if the phases are kinematically unconstrained.
However, for FDMNs with incompressible phases, the online evaluation procedure
needs to be adapted, as we discuss in the following. We define the linear
operator $\underline{\underline{A}}:(\mathds{R}^{3})^{\otimes
N_{\sf{n}}}\rightarrow\left({\sf{Sym}_{0}(3)}\right)^{N_{\sf{w}}}$ which maps
the emulsion-wise jumps $\underline{b}\in(\mathds{R}^{3})^{\otimes
N_{\sf{n}}}$ of the velocity gradient onto the phase-wise strain rate
fluctuations. Also, we express the effective dissipation potential density
$\bar{\Psi}:\left({\sf{Sym}_{0}(3)}\right)^{N_{\sf{w}}}\rightarrow\mathds{R}$
of the FDMN in terms of the phase-wise dissipation potentials
$\Psi_{i}:{\sf{Sym}_{0}(3)}\rightarrow\mathds{R}$, the respective weights
$w_{K+1}^{i}$, and the phase-wise strain rate tensor
$\mbox{\boldmath$D$}_{i}\in{\sf{Sym}_{0}(3)}$ as

$\bar{\Psi}\left(\underline{\bar{D}}+\underline{\underline{A}}\,\underline{b}\right)=\sum_{i=1}^{N_{\sf{w}}}w_{K+1}^{i}\Psi_{i}(\mbox{\boldmath$D$}_{i}),$
(3.27)

where all entries of the vector
$\underline{\bar{D}}\in\left({\sf{Sym}_{0}(3)}\right)^{N_{\sf{w}}}$ equal the
effective strain rate $\bar{\mbox{\boldmath$D$}}$, i.e.,

$\underline{\bar{D}}=(\bar{\mbox{\boldmath$D$}},\bar{\mbox{\boldmath$D$}},...,\bar{\mbox{\boldmath$D$}}).$
(3.28)

Thus, the minimization problem for the online, non-linear evaluation of the
FDMN is

$\bar{\Psi}\left(\underline{\bar{D}}+\underline{\underline{A}}\,\underline{b}\right)\rightarrow\min_{\underline{b}\in(\mathds{R}^{3})^{\otimes
N_{\sf{n}}}},$ (3.29)

and the Euler–Lagrange equation of the problem (3.29) reads

$\partialderivative{\bar{\Psi}}{\underline{b}}\left(\underline{\bar{D}}+\underline{\underline{A}}\,\underline{b}\right)={\underline{\underline{A}}}^{\sf{T}}\partialderivative{\bar{\Psi}}{\underline{D}}\left(\underline{\bar{D}}+\underline{\underline{A}}\,\underline{b}\right)=\underline{0}.$
(3.30)

In terms of a weight matrix
$\underline{\underline{W}}:\left({\sf{Sym}_{0}(3)}\right)^{N_{\sf{w}}}\rightarrow\left({\sf{Sym}_{0}(3)}\right)^{N_{\sf{w}}}$
defined by its action on a vector
${\underline{d}\in\left({\sf{Sym}_{0}(3)}\right)^{N_{\sf{w}}}}$

$\underline{\underline{W}}\,\underline{d}=\left(w_{K+1}^{1}d_{1},w_{K+1}^{2}d_{2},...,w_{K+1}^{N_{\sf{w}}}d_{K}\right),$
(3.31)

we may write the Euler–Lagrange equation (3.30) as

${\underline{\underline{A}}}^{\sf{T}}\underline{\underline{W}}\,\underline{\tau}\left(\underline{\bar{D}}+\underline{\underline{A}}\,\underline{b}\right)=\underline{0},$
(3.32)

with the vector of phase wise stresses
$\underline{\tau}\left(\underline{\bar{D}}+\underline{\underline{A}}\,\underline{b}\right)\in({\sf{Sym}_{0}(3)})^{N}_{\sf{w}}$.
Therefore, the jump vector $\underline{b}$ satisfying equation (3.32) is
determined using a Newton scheme with the update rule

$\underline{b}_{n+1}=\underline{b}_{n}+s_{n}\Delta\underline{b}_{n},$ (3.33)

containing the iteration count $n$, the backtracking factor $s_{n}$, and the
update direction $\Delta\underline{b}_{n}$. The update direction
$\Delta\underline{b}_{n}$ solves the linear system

$\underline{\underline{H}}\Delta\underline{b}_{n}=\left({\underline{\underline{A}}}^{\sf{T}}\underline{\underline{W}}\,\underline{\tau}\left(\underline{\bar{D}}+\underline{\underline{A}}\,\underline{b}\right)\underline{\underline{A}}\right)\Delta\underline{b}_{n}={\underline{\underline{A}}}^{\sf{T}}\underline{\underline{W}}\,\underline{\tau}\left(\underline{\bar{D}}+\underline{\underline{A}}\,\underline{b}\right).$
(3.34)

So far, the procedure is completely analogous to the approaches for (direct)
DMNs presented in previous work [54, 46, 59]. However, in the case of
incompressible phases, the matrix $\underline{\underline{A}}$ is symmetric but
singular, and the solution of the equation system $\eqref{eq:updaterule}$ is
not unique. However, because of incompressibility, we know that the $i$-th
entry $\mbox{\boldmath$b$}_{i}$ of the jump vector $\underline{b}$ is
orthogonal to the $i$-th entry $\mbox{\boldmath$n$}_{i}$ of the vector of
layering normals $\underline{n}$, i.e.,

$\mbox{\boldmath$b$}_{i}\cdot\mbox{\boldmath$n$}_{i}=0.$ (3.35)

Hence, we construct the block diagonal matrix
$\underline{\underline{N}}:(\mathds{R}^{3})^{\otimes
N_{\sf{n}}}\rightarrow(\mathds{R}^{3})^{\otimes N_{\sf{n}}}$ consisting of
$N_{\sf{w}}$ blocks $\mbox{\boldmath$N$}_{i}$ on the main diagonal, which are
defined as

$\mbox{\boldmath$N$}_{i}=\beta_{i}\mbox{\boldmath$n$}_{i}\otimes\mbox{\boldmath$n$}_{i},\quad
i\in\\{1,...,N_{\sf{w}}\\}$ (3.36)

with the constant $\beta_{i}\in\mathds{R}_{>0}$. By definition it holds that

$\mbox{\boldmath$N$}_{i}\mbox{\boldmath$b$}_{i}=\mbox{\boldmath$0$},$ (3.37)

and the jump vector $\underline{b}$ is in the kernel of
$\underline{\underline{N}}$, such that

$\underline{\underline{N}}\,\underline{b}=\underline{0}.$ (3.38)

Thus, we add equations (3.34) and (3.38), such that the solution
$\Delta\underline{b}_{n}$ of the resulting update rule

$\left(\underline{\underline{H}}+\underline{\underline{N}}\right)\Delta\underline{b}_{n}={\underline{\underline{A}}}^{\sf{T}}\underline{\underline{W}}\,\underline{\tau}\left(\underline{\bar{D}}+\underline{\underline{A}}\,\underline{b}\right),$
(3.39)

is unique. However, floating point precision may cause issues during the
addition of the operators $\underline{\underline{H}}$ and
$\underline{\underline{N}}$, if numbers of largely different magnitudes are
summed. To remedy this issue, the scalars $\beta_{i}$ need to be chosen
according to the magnitude of the entries in the matrix
$\underline{\underline{H}}$. Furthermore, some input weights $w_{K+1}^{i}$
might become equal to, or close to zero during training, which might render
the system matrix $(\underline{\underline{H}}+\underline{\underline{N}})$ ill-
conditioned. Thus, it might be difficult to obtain accurate update directions
$\Delta\underline{b}_{n}$ by solving equation (3.39). While it is generally
possible to prune the FDMN tree [53] by removing phases with vanishing or
almost vanishing volume fractions, this might collapse the rank-$3$ CLMs of
the bottom FDMN layer into layered materials of a lesser rank. In case of
infinite material contrast, this could result in the propagation of singular
effective properties through the FDMN. It is therefore preferable to work with
possibly ill-conditioned systems and to employ appropriate methods for the
solution of ill-conditioned systems [78]. Alternatively, multiple FDMNs can be
trained and the ones with the best online evaluation results can be selected.
In this article, we follow the latter approach.

##  4 Application to rigid fibers suspended in polyamide 6

###  4.1 Material description and computational aspects

Let a fiber orientation state be defined by the fiber orientation distribution
function $\rho:\mathcal{S}^{2}\rightarrow\mathds{R}$, which encodes the
probability that fibers are oriented in the direction
$\mbox{\boldmath$p$}\in\mathcal{S}^{2}$. In component scale molding
simulations, the orientation distribution function $\rho$ is often spatially
inhomogeneous, and it is therefore prohibitively expensive to compute the
spatial and temporal evolution of the function $\rho$ [79, 80]. Therefore, the
evolution of surrogate quantities that contain some of the information encoded
in the function $\rho$ is often considered instead. One such quantity is the
second order fiber orientation tensor [81, 82]

$\mbox{\boldmath$N$}=\int_{\mathcal{S}^{2}}\mbox{\boldmath$p$}\otimes\mbox{\boldmath$p$}\,\rho(\mbox{\boldmath$n$})\mathop{}\\!\mathrm{d}S(\mbox{\boldmath$n$}),$
(4.1)

which is commonly used in component scale molding simulations to track the
evolution of the fiber orientation state [79, 80].

(a) Fiber orientation triangle $S_{\sf{T}}$

(b) Cross-type fit

Figure 5: Fiber orientation triangle $S_{\sf{T}}$ in CMYK coloring with 31
evaluation points (a), and material data with Cross-type fit for Ultramid®B3K
(b).

Because a major prospect of FDMNs is the combination of component scale
molding simulations with high fidelity information from full field viscous
computations, we employ the orientation tensor $N$ to describe the fiber
orientation state of the considered microstructures. The tensor $N$ may be
parametrized by its two largest eigenvalues $\lambda_{1}$ and $\lambda_{2}$,
as well as a proper orthogonal tensor $R$ via [83, 84]

$\mbox{\boldmath$N$}=\mbox{\boldmath$R$}\,\mbox{\boldmath$\mathrm{diag}$}(\lambda_{1},\lambda_{2},1-\lambda_{1}-\lambda_{2})\,{\mbox{\boldmath$R$}}^{\sf{T}},\quad\lambda_{1}\geq\lambda_{2}.$
(4.2)

By objectivity [51], and because the tensor $N$ has unit trace and is positive
semi-definite, every fiber orientation tensor $N$ represents a point in the
fiber orientation triangle $S_{\sf{T}}$ [55, 84]

$S_{\sf{T}}=\left\\{\underline{\lambda}={(\lambda_{1},\lambda_{2})}^{\sf{T}}\quad\middle|\quad
1\geq\lambda_{1}\geq\frac{1}{3}\quad\textrm{and}\quad\min\left(\lambda_{1},1-\lambda_{1}\right)\geq\lambda_{2}\geq\frac{1-\lambda_{1}}{2}\right\\}.$
(4.3)

For detailed discussions regarding the description of fiber orientation
states, we refer the reader to Kanatani [81] and Bauer et al. [84, 85]. To
generate the training data and the online validation data for the FDMN, we use
FFT-based computational procedures to compute the effective viscous response
of fiber suspensions [38, 17]. We generate the required fiber suspension
microstructures using the sequential addition and migration method [86] for 31
orientation states $\underline{\lambda}$ in the fiber orientation triangle
$S_{\sf{T}}$, see Figure 7. For the visualization of the orientation triangle
$S_{\sf{T}}$, we follow Köbler et al. [55], and use a CMYK coloring scheme. To
keep the computational effort for computing the effective viscous responses
feasible, we restrict to microstructures with a fiber volume fraction of 25%,
and assign all fibers the length $\ell$ and the diameter $d$. Also, we set the
aspect ratio $r_{a}=\ell/d$ of all fibers to 10, and use the results of Sterr
et al. [17] regarding the required resolution and size of the microstructure
volume elements. Consequently, we choose a resolution of 15 voxels per fiber
diameter, as well as cubic volume elements with edge size $L=2.2\ell$.
Regarding the matrix material, we build on the investigations in Sterr et al.
[17, 51], and prescribe a commercially available polyamide 6 [87] as the
matrix material. We fit a Cross-type constitutive equation

$\eta(\dot{\gamma})=\eta_{\infty}+\frac{\eta_{0}-\eta_{\infty}}{1+(k\dot{\gamma})^{m}},$
(4.4)

to the available material data in the interval
${[1.7,16300]\,\mathrm{s}^{-1}}$ at a temperature of
${250^{\circ}\mathrm{C}}$, and collect the resulting model parameters in Table
1.

$\eta_{0}$ | $\eta_{\infty}$ | $k$ | $m$
---|---|---|---
288.9 $\mathrm{Pa}\,\mathrm{s}$ | 15.0 $\mathrm{Pa}\,\mathrm{s}$ |  $10.9\cdot 10^{-4}$ $\mathrm{s}$ | 1.1
Table 1: Parameters of the Cross-type constitutive equation (4.4) for a
commercially available polyamide 6 [87].

For all non-Newtonian considerations, we investigate the six load cases
collected in the matrix $\underline{\underline{\bar{D}}}$ in Mandel notation

$\underline{\underline{\bar{D}}}=\dot{\gamma}\sqrt{\frac{2}{3}}\left[\begin{array}[]{rrrrrr}1&-\frac{1}{2}&-\frac{1}{2}&0&0&0\\\
-\frac{1}{2}&1&-\frac{1}{2}&0&0&0\\\ -\frac{1}{2}&-\frac{1}{2}&1&0&0&0\\\
0&0&0&\sqrt{\frac{3}{2}}&0&0\\\ 0&0&0&0&\sqrt{\frac{3}{2}}&0\\\
0&0&0&0&0&\sqrt{\frac{3}{2}}\end{array}\right],$ (4.5)

or each macroscopic scalar shear rate $\dot{\gamma}$ in the set of studied
shear rates $S_{\dot{\gamma}}$, such that

$\dot{\gamma}\in S_{\dot{\gamma}}={\left\\{a\cdot 10^{b}\ \mathrm{s^{-1}}|\
a=1,2,5;\ b=1,2,3,4\right\\}\cup\left\\{10^{5}\,\mathrm{s}^{-1}\right\\}}.$
(4.6)

Consequently, the set $\bar{\mbox{\boldmath$D$}}_{\dot{\gamma}}$ of all
investigated load cases is defined by the equations (4.5) and (4.6). Because
the viscous stress inside the rigid fibers is, constitutively, not well-
defined [17], we employ a dual formulation of the associated homogenization
problem for the FFT-based computations. We discretize the microstructures on a
staggered grid [88], and solve the resulting equation systems with the
conjugate gradient method (CG) for linear computations, and with a Newton-CG
approach for non-linear computations. With the goal of training an FDMN for
each of the 31 fiber orientation states shown in Figure 5(a), we generate 32
sample viscosities with the procedure described in section 3.3, and compute
the corresponding effective viscosities. To compute a single effective
viscosity, five FFT-based computations are necessary [38], leading to a total
of 4960 computations, and 992 computed effective viscosities. In the previous
sections 2.3 and 3.2, we studied the linear homogenization functions of CLMs,
and presented the FDMN architecture to treat infinite material contrast and
incompressible materials. In the following, we discuss the material sampling
for the offline training, the offline training procedure, and the online
evaluation of FDMNs for suspensions of rigid particles.

###  4.2 Offline training

With the linear training data at hand, we wish to train FDMNs for the
prediction of the non-Newtonian viscous behavior of shear-thinning fiber
polymer suspensions. We follow previous work [54, 46] and choose the depth of
the FDMN as $K=8$ to achieve sufficient prediction quality for non-linear
computations. An FDMN of depth $K=8$ with CLMs of rank $R=3$ has 512 weights
and 638 angles as free parameters, as per equations (3.6) and (3.20). For each
microstructure, we define the training data $\mathcal{S}^{\sf{D}}$ as

$\mathcal{S}^{\sf{D}}=\\{({\mathbb{V}}^{s}_{1},{\mathbb{V}}_{F},\bar{{\mathbb{V}}}^{s})\quad|\quad
s\in(1,...,32)\\},$ (4.7)

where ${\mathbb{V}}^{s}_{1}$ denotes the sample viscosities generated with
equation (3.7), ${\mathbb{V}}_{F}$ denotes the infinite viscosity of the rigid
fibers, and $\bar{{\mathbb{V}}}^{s}$ stands for the effective viscosity of the
sample $s$. For each microstructure, we split the training data
$\mathcal{S}^{\sf{D}}$ into the training set
$\mathcal{S}^{\sf{t}}\subset\mathcal{S}^{\sf{D}}$ and the validation set
$\mathcal{S}^{\sf{v}}\subset\mathcal{S}^{\sf{D}}$, which consist of 90% and
10% of the total training data, respectively. The training set
$\mathcal{S}^{\sf{t}}$ and the validation set $\mathcal{S}^{\sf{v}}$ share no
samples, i.e., $\mathcal{S}^{\sf{t}}\cap\mathcal{S}^{\sf{v}}=\emptyset$.

As an initial guess for the FDMN parameters, we uniformly sample all angles
$\underline{\alpha}$ and weights $\underline{w}$ from their respective
intervals, and rescale the weights $\underline{w}$ such that they sum to
unity. Because we use a relatively small set of 32 samples as training data,
we dedicate a large portion of it to the training set $\mathcal{S}^{\sf{t}}$,
and choose to train 20 FDMNs per fiber orientation state
$\underline{\lambda}$. Furthermore, we train the FDMNs on mini batches with
size $N_{\sf{b}}=8$, for which we draw randomly from the training set
$\mathcal{S}^{\sf{t}}$. In case the last batch is smaller than eight samples,
we drop the batch. Because the initialization of the parameter vectors
$\underline{\alpha}$ and $\underline{w}$, as well as the offline training
process are random, the FDMNs differ in their parameters and quality of fit.
By training 20 FDMNs per fiber orientation state, we leverage this randomness
with the goal of obtaining FDMNs with a high quality of fit. This strategy
aims to reduce the total computational effort required to obtain sufficiently
accurate FDMNs, because repeatedly training FDMNs on 32 samples requires less
computational resources than conducting full field simulations for larger
sample sizes. For the minimization of the loss function $\mathcal{L}$, we
choose the RAdam algorithm [77] and use a learning rate sweep [89] to
determine the learning rates $\kappa_{0,\alpha}$ and $\kappa_{0,w}$. With the
learning rate sweep we obtain highly similar learning rates for both parameter
groups $\underline{\alpha}$ and $\underline{w}$, such that
${\kappa_{0,\alpha}=\kappa_{0,w}=1\cdot 10^{-2}}$. As the learning process
advances, we use PyTorch’s StepLR learning rate scheduler to improve
convergence towards minima [90, 91]. To do so, we multiply the learning rates
$\kappa_{0,\alpha}$ and $\kappa_{0,w}$ with a constant factor ${f_{LR}=0.75}$
every 150 epochs. Other than the learning rates, we used standard
hyperparameters for the two momentum coefficients $\beta_{1}=0.9$ and
$\beta=0.999$, as well as the stabilization constant $\varepsilon=10^{-8}$.
Overall, we train every FDMNs for a total of 2000 epochs each.

(a) Training set $\mathcal{S}^{\sf{t}}$

(b) Validation set $\mathcal{S}^{\sf{v}}$

Figure 6: Minimum, maximum, and average mean error $e^{\sf{mean}}$ of all
trained FDMNs for the training set $\mathcal{S}^{\sf{t}}$ (a) and the
validation set $\mathcal{S}^{\sf{v}}$ (b).

To measure the quality of fit, we define the mean error $e^{\sf{mean}}$ over a
subset $\mathcal{S}$ of the training data $\mathcal{S}^{\sf{D}}$ as

$e^{\sf{mean}}=\frac{1}{|\mathcal{S}|}\sum_{\left({\mathbb{V}}_{1}^{s},{\mathbb{V}}_{F},\bar{{\mathbb{V}}}^{s}\right)\in\,\mathcal{S}}\frac{\norm*{\mathcal{DMN}_{\sf{L}}({\mathbb{V}}_{1}^{s},{\mathbb{V}}_{F},\underline{\alpha},\underline{w})-\bar{{\mathbb{V}}}^{s}}_{1}}{\norm*{\bar{{\mathbb{V}}}^{s}}_{1}},$
(4.8)

where $|\mathcal{S}|$ is the cardinality of the set $\mathcal{S}$. Depending
on the fiber orientation state $\underline{\lambda}$ and the initial guesses
of the angles $\underline{\alpha}$ and the weights $\underline{w}$, the
quality of fit varies per DMN. Therefore, for the whole training process of
all 620 FDMNs, we visualize the smallest mean error $\min(e^{\sf{mean}})$, the
largest mean error $\max(e^{\sf{mean}})$, and the average mean error
$\mathrm{avg}(e^{\sf{mean}})$ for the training and validation sets
$\mathcal{S}^{\sf{t}}$ and $\mathcal{S}^{\sf{v}}$ per epoch in Figure 6. For
the training and the validation set, we observe that the three considered
error measures drop rapidly at the beginning of the training, and continue to
improve as the training continues. As the learning rates are small for epochs
larger than 1500, the errors do not change significantly at the end of the
training, and convergence is ensured. Evidently, the largest mean error
$\max(e^{\sf{mean}}_{\sf{t}})$ on the training set $\mathcal{S}^{\sf{t}}$
drops below 1% at the end of the training process, while the smallest mean
error $\min(e^{\sf{mean}}_{\sf{t}})$ falls below 0.3%. This indicates a high
quality of fit on the training set $\mathcal{S}^{\sf{t}}$ for all FDMNs.
Compared to the training set $\mathcal{S}^{\sf{t}}$, the spread between the
largest mean error $\max(e^{\sf{mean}}_{\sf{v}})$ and the smallest mean error
$\min(e^{\sf{mean}}_{\sf{v}})$ for the validation set $\mathcal{S}^{\sf{v}}$
is larger. However, the largest mean error $\max(e^{\sf{mean}}_{\sf{v}})$
stays below 2% for all considered FDMNs and does not fluctuate by a large
margin. Additionally, the average mean error
$\mathrm{avg}(e^{\sf{mean}}_{\sf{v}})$ decreases continuously during training.
Consequently, the prediction quality on both the training and validation set
improve on average as the training progresses, although the largest mean
validation error $\max(e^{\sf{mean}}_{\sf{v}})$ does not improve significantly
after 1250 training epochs. Overall, this indicates that no pronounced
overfitting to the training data set $\mathcal{S}^{\sf{t}}$ occurs.

###  4.3 Online evaluation

To measure the online performance of an FDMN, we define the online error
function $e_{\sf{on}}:{\sf{Sym}_{0}(3)}\rightarrow\mathds{R}$ as

$e_{\sf{on}}(\bar{\mbox{\boldmath$D$}})=\frac{\norm*{\mbox{\boldmath$\sigma$}^{DMN}(\bar{\mbox{\boldmath$D$}})-\mbox{\boldmath$\sigma$}^{FFT}(\bar{\mbox{\boldmath$D$}})}_{2}}{\norm*{\mbox{\boldmath$\sigma$}^{FFT}(\bar{\mbox{\boldmath$D$}})}_{2}},$
(4.9)

where $\norm*{\cdot}_{2}$ denotes the $\ell^{2}$-norm of the components in
Mandel notation, $\bar{\mbox{\boldmath$D$}}$ stands for the prescribed
effective strain rate tensor, and
$\mbox{\boldmath$\sigma$}^{DMN}:{\sf{Sym}_{0}(3)}\rightarrow{\sf{Sym}_{0}(3)}$
and
$\mbox{\boldmath$\sigma$}^{FFT}:{\sf{Sym}_{0}(3)}\rightarrow{\sf{Sym}_{0}(3)}$
are the stress functions of the FDMN and the FFT-based homogenization,
respectively. We then choose the largest error over all load cases
$e_{\sf{on}}^{\sf{max}}$ and the mean error over all load cases
$e_{\sf{on}}^{\sf{mean}}$

$e_{\sf{on}}^{\sf{max}}=\max_{\bm{\bar{D}}\in\bm{\bar{D}}_{\dot{\gamma}}}e_{\sf{on}},\quad
e_{\sf{on}}^{\sf{mean}}=\frac{1}{|\bar{\mbox{\boldmath$D$}}_{\dot{\gamma}}|}\sum_{\bm{\bar{D}}\in\bm{\bar{D}}_{\dot{\gamma}}}e_{\sf{on}}(\bar{\mbox{\boldmath$D$}}),$
(4.10)

to study the performance of the different FDMNs. For each investigated
orientation state $\underline{\lambda}$, we identify the FDMN with the
smallest maximum error $e_{\sf{on}}^{\sf{max}}$ and visualize the associated
errors $e_{\sf{on}}^{\sf{mean}}$ and $e_{\sf{on}}^{\sf{max}}$ over the fiber
orientation triangle, see Figures 7(a) and 7(b)

(a) Mean error $e_{\sf{on}}^{\sf{mean}}$

(b) Largest error $e_{\sf{on}}^{\sf{max}}$

Figure 7: Mean online error $e_{\sf{on}}^{\sf{mean}}$ (a) and largest online
error $e_{\sf{on}}^{\sf{max}}$ (b) for the best FDMNs over the fiber
orientation triangle.

The maximum error $e_{\sf{on}}^{\sf{max}}$ of the best FDMNs range between
2.88% for the orientation state
${\underline{\lambda}={(0.4444,0.3611)}^{\sf{T}}}$, and 4.31% for the fiber
orientation state $\underline{\lambda}={(0.8055,0.1389)}^{\sf{T}}$. Thus, for
the investigated orientation states, the FDMNs yield sufficient prediction
accuracy in terms of engineering requirements over a wide range of shear rates
$\dot{\gamma}\in[10,10^{5}]\,\mathrm{s}^{-1}$. Additionally, the relatively
low maximum errors $e_{\sf{on}}^{\sf{max}}$ show that with fewer samples than
in previous studies for solid materials, for example by Liu et al. [52, 53] or
Gajek et al. [54, 46, 59], appropriately accurate FDMNs can be produced for
suspensions of rigid fibers. Generally, the maximum error
$e_{\sf{on}}^{\sf{max}}$ is larger for more strongly aligned oriented
orientation states towards the lower right-hand side of the fiber orientation
triangle than for less strongly aligned oriented orientation states towards
the upper and the left-hand side. The mean error $e_{\sf{on}}^{\sf{mean}}$
ranges between 1.16% for the orientation state
$\underline{\lambda}={(0.8334,0.0833)}^{\sf{T}}$, and 1.82% for the
orientation state $\underline{\lambda}={(0.6528,0.3194)}^{\sf{T}}$,
underlining the prediction accuracy of the FDMNs.

To discuss one specific example, we take a more detailed look at the
performance of the best FDMN at the fiber orientation state
$\underline{\lambda}={(0.8055,0.1389)}^{\sf{T}}$ where the largest maximum
error $e_{\sf{on}}^{\sf{max}}=4.31\%$ occurs. Let the components of the strain
rate tensor $\tilde{\mbox{\boldmath$D$}}_{i}\in{\sf{Sym}_{0}(3)}$, in the
standard basis $\\{\mbox{\boldmath$e$}_{i}\\}$ and in Mandel notation, be
given by the $i$-th column of the matrix
$\underline{\underline{\bar{D}}}/\dot{\gamma}$, see equation (4.5). Then, the
tensor $\tilde{\mbox{\boldmath$D$}}_{i}$ represents a load with unit norm, and
we refer to the tensor $\tilde{\mbox{\boldmath$D$}}_{i}$ as load case. In the
following, we discuss the observed stress responses and the prediction quality
of the FDMN for all considered shear rates (4.6) and the six load cases
defined in equation (4.5), see Figure 8. Because of the strong alignment of
the fibers in the coordinate direction $\mbox{\boldmath$e$}_{1}$, the stress
norms are the largest for the load case $\tilde{\mbox{\boldmath$D$}}_{1}$,
which encodes incompressible elongational flow in the coordinate direction
$\mbox{\boldmath$e$}_{1}$, see Figure 8(a).

(a) $\ell^{2}$-norm of the effective stress response
$\bar{\mbox{\boldmath$\sigma$}}(\bar{\mbox{\boldmath$D$}})$

(b) Online error $e_{\sf{on}}(\bar{\mbox{\boldmath$D$}})$

Figure 8: $\ell^{2}$-norm of the effective stress response
$\bar{\mbox{\boldmath$\sigma$}}(\bar{\mbox{\boldmath$D$}})$ over the shear
rate $\dot{\gamma}\in[10,10^{5}]\,\mathrm{s}^{-1}$ at the orientation state
${\underline{\lambda}=(0.8055,0.1389)}$ as computed with FFT-based
homogenization and the best identified FDMN (a). Online error
$e_{\sf{on}}(\bar{\mbox{\boldmath$D$}})$ of the FDMN over the shear rate
$\dot{\gamma}\in[10,10^{5}]\,\mathrm{s}^{-1}$ (b).

The norms of the stress responses to elongational flow in the coordinate
directions $\mbox{\boldmath$e$}_{2}$ and $\mbox{\boldmath$e$}_{3}$, i.e., the
load cases $\tilde{\mbox{\boldmath$D$}}_{2}$ and
$\tilde{\mbox{\boldmath$D$}}_{3}$, are comparatively small. For the largest
considered shear rate of $\dot{\gamma}=10^{5}\,\mathrm{s}^{-1}$, we observe a
stress norm of 21 $\mathrm{MPa}$ for the load case
$\tilde{\mbox{\boldmath$D$}}_{1}$, as well as 11.8 $\mathrm{MPa}$ and 12.1
$\mathrm{MPa}$ for the load cases $\tilde{\mbox{\boldmath$D$}}_{2}$ and
$\tilde{\mbox{\boldmath$D$}}_{3}$, respectively. Because of the similar
degrees of fiber alignment in the coordinate directions
$\mbox{\boldmath$e$}_{2}$ and $\mbox{\boldmath$e$}_{3}$, the stress norms
associated with the load cases $\tilde{\mbox{\boldmath$D$}}_{2}$ and
$\tilde{\mbox{\boldmath$D$}}_{3}$ differ only slightly. Furthermore, the
stress norms in response to shear in the load cases
$\tilde{\mbox{\boldmath$D$}}_{4}$ and $\tilde{\mbox{\boldmath$D$}}_{5}$ are
smaller than to shear in the load cases $\tilde{\mbox{\boldmath$D$}}_{6}$.
Again, this is caused by the strong alignment of the fibers in the coordinate
direction $\mbox{\boldmath$e$}_{1}$, which increases the flow resistance
stronger in the $\mbox{\boldmath$e$}_{1}$-$\mbox{\boldmath$e$}_{2}$-plane than
in other shear planes [51]. For the largest considered shear rate of
$\dot{\gamma}=10^{5}\,\mathrm{s}^{-1}$, the stress norms are 5.5
$\mathrm{MPa}$, 6.0 $\mathrm{MPa}$, and 8.4 $\mathrm{MPa}$ for the load cases
$\tilde{\mbox{\boldmath$D$}}_{4}$, $\tilde{\mbox{\boldmath$D$}}_{5}$, and
$\tilde{\mbox{\boldmath$D$}}_{6}$, respectively. Overall, the predicted
$\ell^{2}$-norm of the effective stress response
$\bar{\mbox{\boldmath$\sigma$}}^{DMN}$ of the FDMN agrees closely with the
effective stress $\bar{\mbox{\boldmath$\sigma$}}^{FFT}$ computed via FFT-based
homogenization for all investigated load directions
$\tilde{\mbox{\boldmath$D$}}_{i},\,i\in\\{1,...,6\\}$, and shear rates
$\dot{\gamma}\in S_{\dot{\gamma}}$, see Figure 8(a).

Not only the $\ell^{2}$-norm of the effective stress response is predicted
well by the FDMN, but the direction as well. This is evident from the
relatively small online errors $e_{\sf{on}}(\bar{\mbox{\boldmath$D$}})$ shown
in 8(b). Depending on the load case, the error $e_{\sf{on}}$ increases up to a
point in the interval $[10^{2},10^{3}]\,\mathrm{s}^{-1}$, before decreasing
again for higher shear rates. The largest error occurs for the elongational
load case $\tilde{\mbox{\boldmath$D$}}_{3}$ with 4.31%, and is closely
followed by 4.08% for the elongational load case
$\tilde{\mbox{\boldmath$D$}}_{1}$. For the shear load cases
$\tilde{\mbox{\boldmath$D$}}_{4}$, $\tilde{\mbox{\boldmath$D$}}_{5}$, and
$\tilde{\mbox{\boldmath$D$}}_{6}$, the largest errors are 3.05%, 2.78%, and
3.23%, respectively, while the largest error for the elongational load case
$\tilde{\mbox{\boldmath$D$}}_{2}$ is 2.2%. The largest error for the shear
load case $\tilde{\mbox{\boldmath$D$}}_{6}$ occurs at a shear rate of
$\dot{\gamma}=10^{2}\,\mathrm{s}^{-1}$, whereas the largest errors for all
other load cases occur at a shear rate of $\dot{\gamma}=5\cdot
10^{2}\,\mathrm{s}^{-1}$. This is directly tied to the Cross-type constitutive
equation (4.4), which has two Newtonian plateaus for low and high shear rates
and a non-linear transition in between. Consequently, for large shear rates
$\dot{\gamma}\rightarrow\infty$ the Matrix behavior is Newtonian and the
homogenization function of the microstructure approaches the linear
homogenization function. Since the FDMN approximates the non-linear
homogenization function the microstructure to first order [54], the prediction
quality depends on the degree to which the matrix behaves non-linearly. This
is in line with previous observations regarding the modelling of shear-
thinning fiber suspensions [17, 51].

###  4.4 Computational Speedup

| Sampling total | Sampling per orientation (mean) | Training
---|---|---|---
CPU Hours | 24415 | 788 | 0.08
Wall-clock time (h) | 305 | 9.8 | 0.08
Table 2: CPU hours and wall-clock time for the training of one FDMN and the sampling of effective viscosities $\bar{{\mathbb{V}}}^{s}$ using FFT-based homogenization. | FFT (h) | DMN (ms) | Speedup
---|---|---|---
Min | 2.08 | 600 | 11785
Mean | 2.61 | 630 | 14828
Max | 3.07 | 650 | 17225
Table 3: Speedups and wall-clock times for the evaluation of a single load
case using FFT-based homogenization compared to an FDMN.

For sampling the linear homogenization functions of the 31 considered fiber
orientation states $\underline{\lambda}\in S_{\sf{T}}$, we used a workstation
with two AMD EPYC 7552 48-Core processors and 1024 GB DRAM. We relied on a
single thread per FFT-based computation and ran 80 computations in parallel.
In total, the computations took 24415 CPU hours and a wall-clock time of 305
hours, see Table 2. Averaging over all 31 considered orientations
$\underline{\lambda}$, this leads to 788 CPU Hours and $9.8\mathrm{h}$ wall-
clock time per orientation. For the training and the online evaluation of the
FDMNs, as well as the non-linear FFT-based computations, we used a workstation
with two AMD EPYC 9534 64-Core Processors and 1024 GB DRAM, and ran all
computations on a single thread. The wall-clock time to evaluate all six load
cases $\tilde{\mbox{\boldmath$D$}}_{i},\,i=1,...,6$, and thirteen shear rates
$\dot{\gamma}\in S_{\dot{\gamma}}$ in series for a single fiber orientation
state $\underline{\lambda}$ is listed in Table 3. The wall-clock times for the
FFT-based computations range from $2.08\mathrm{h}$ to $3.07\mathrm{h}$ with
the mean wall-clock time over all considered orientation states
$\underline{\lambda}$ being $2.61\mathrm{h}$. In contrast, the online
evaluation of the FDMNs takes between $600\mathrm{ms}$ and $630\mathrm{ms}$,
with a mean of $650\mathrm{ms}$. Consequently, the speedup factors range
between 11785 and 17225 with a mean of 14828. This considerable speedup is
achieved by investing computational resources into the sampling of the linear
homogenization functions of the microstructures. It is straightforward to
judge whether this investment is sensible for the considered setup. For one
microstructure the generation of training data and the training of a single
FDMN takes $9.88\mathrm{h}$ wall-clock time on average, while the non-linear
FFT-based computations take $2.61\mathrm{h}$ wall-clock time on average.
Therefore, assuming the non-linear computations run on one thread as in the
presented setup, it would take four non-linear computation to offset the
initial investment into an FDMN on average. Given that component scale
simulations routinely require thousands or millions of microscale computations
[46, 59], the sampling and training effort is offset easily in engineering
problems.

###  4.5 Comparison with machine learning aided analytical models

In this section, we compare the FDMN based approach to predicting the
effective behavior of fiber suspensions with another machine learning approach
suggested by Sterr et al. [51]. The authors presented four different
analytical constitutive equations for the effective viscosity of shear-
thinning fiber suspensions, and identified the model parameters using
supervised machine learning. The model parameters were learned from simulation
data obtained with FFT-based computational techniques for the same material as
considered in this article: rigid fibers suspended in a Cross-type matrix
material with parameters as shown in Table 1. However, the FDMN approach and
the analytical approach follow very different paradigms. The FDMNs are trained
to approximate the non-linear homogenization function of the microstructure by
learning the respective linear homogenization function without prior knowledge
of the actual constitutive equations. In contrast, the parameters of the
analytical models were learned from the non-linear stress response of the
suspension, and the analytical models incorporate knowledge about the expected
material behavior. Even though the two approaches follow different strategies,
they achieve comparable prediction accuracy. For the load cases
$\bar{\mbox{\boldmath$D$}}_{\dot{\gamma}}$ and fiber orientation states
$\underline{\lambda}$ considered in this article, the largest validation error
of the FDMNs is $4.31\%$ and occurs for the fiber orientation state
$\underline{\lambda}={(0.8055,0.1389)}^{\sf{T}}$. For the same load cases
$\bar{\mbox{\boldmath$D$}}_{\dot{\gamma}}$ and the same set of fiber
orientation states $\underline{\lambda}$, three of the four analytical models
by Sterr et al. [51] achieve a maximum prediction error of $5.00\%$, which
occurs for the fiber orientation state
${\underline{\lambda}={(0.8334,0.0833)}^{\sf{T}}}$. The fourth model did not
compare favorably with the other three, because the built-in assumptions of
stress-strain rate superposition and orientation averaging did not hold well
for the considered suspension [51]. Consequently, the FDMNs achieved a
slightly higher prediction accuracy for the type of fiber suspension
considered in this article. Also, within the constraints of a first order
approximation, FDMNs generalize to different constitutive equations while the
presented analytical models do not, reducing the required modelling effort for
FDMNs. However, the FDMNs utilized in this article use 1150 free parameters,
whereas the analytical models use between 11 and 49 parameters.

In terms of computational cost, generating the training data and training an
FDMN for a single microstructure took $9.88\mathrm{h}$ wall-clock time on
average, see Table 2. The parameters of the analytical models for a single
microstructure were identified based on six non-linear FFT-based simulations,
which took $2.61\mathrm{h}$ wall-clock time on average. This leads to an
average total wall-clock time of $15.70\mathrm{h}$ if all required simulations
are run in series on one thread, and $0.04\,\mathrm{h}$ are allocated to
identify the model parameters. However, with the considered setup, the non-
linear FFT-based simulations could be run in parallel, such that the
analytical model parameters can be obtained more quickly than an FDMN. In
summary, FDMNs offer a higher degree of accuracy and flexibility regarding the
considered constitutive equations, but require more computational resources
and a larger amount of parameters. Furthermore, FDMNs are inherently
thermodynamically consistent [54, §3.1], and inherit stress-strain rate
monotonicity from their phases [54, App. C]. In contrast, both thermodynamic
consistency and stress-strain rate monotonicity must be ensured manually in
constitutive models.

##  5 Conclusions

In this article, we extended the direct DMN architecture to the Flexible DMN
(FDMN) architecture for the treatment of fiber suspensions with infinite
material contrast and shear-thinning matrix behavior. To do so, we derived
linear homogenization functions for two-phase layered emulsions that are
governed by Stokes flow and consist of linearly viscous phases. More
specifically, we utilized results by Kabel et al. [70] and Milton [69] on the
homogenization of laminates to derive closed form analytical expressions for
the effective properties of such layered emulsions. Because rank-one layered
materials are ill-suited as DMN building blocks in case of infinite material
contrast, we investigated the effective properties of coated layered materials
(CLMs). We leveraged the linear homogenization functions of rank-one layered
materials and investigated under which conditions the effective behavior of
CLMs is non-singular if the core material is rigid. The conditions depend on
the physical constraints of the employed materials, and involve the rank of
the CLM and the relative orientation of the layering directions. In the
relevant case for incompressible fiber suspensions, a CLM consists of an
incompressible coating material and a rigid core material. Then, the effective
viscosity of a rank-$3$ CLM is non-singular if the three layering directions
are mutually non-orthogonal and mutually non-collinear.

Using the derived homogenization functions for layered materials and CLMs, we
extended the (direct) DMN architecture to the FDMN architecture by replacing
the lowest layer of rank-one laminates in a (direct) DMN with non-singular
CLMs, and the other rank-one laminates with rank-one layered materials capable
of treating fluids. For the offline training of an FDMN, we presented a
strategy where the relative angles of the CLM layering directions are fixed to
reduce the amount of free parameters and guarantee non-singular CLMs.
Furthermore, we modified the online evaluation strategy of direct DMNs to
account for incompressible phases. We leveraged the FDMN architecture to
predict the non-linear effective behavior of fiber suspensions with a Cross-
type matrix material. Compared to direct numerical simulations with FFT-based
computational techniques, the FDMNs achieved validation errors below $4.31\%$
for a variety of 31 fiber orientation states, six different load cases, and a
wide range of shear rates relevant to engineering processes. If the time
required to generate the training data and train the FDMNs is not considered,
the FDMNs achieved an average speed up factor of 14828 as compared to FFT-
based simulations. Additionally, FDMNs achieve higher accuracy than another
machine learning based approach by Sterr et al. [51] for the same composite
material. However, this competitive accuracy improvement comes at the cost of
an increased computational effort to obtain the training data for the FDMNs.

In future work, the application and extensions of the FDMN architecture to
problems involving more than two phases, where multiple phases may be
singular, could be explored. This could enable the FDMN based treatment of
three-phase systems containing a fluid, fibers, and gas, which occur in
engineering processes like injection molding [92] and flotation de-inking
[93]. Furthermore, problems involving solids with defects could be explored as
well. To account for the locally varying microstructure in component scale
simulations, fiber orientation interpolation schemes [57, 46, 63] could be
combined with the FDMN architecture to increase versatility, and enable new
concurrent two-scale simulations involving defects and rigid inclusions.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## Data availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

## Acknowledgements

The research documented in this manuscript was funded by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation), project number
255730231, within the International Research Training Group “Integrated
engineering of continuous-discontinuous long fiber-reinforced polymer
structures“ (GRK 2078/2). The support by the German Research Foundation (DFG)
is gratefully acknowledged.

## Author contributions

The present study was conceptualized by B. Sterr, M. Schneider and T. Böhlke.
The presented FDMN architecture was derived by B. Sterr, S. Gajek, M.
Schneider, and T. Böhlke. B. Sterr, S. Gajek, and M. Schneider implemented and
validated the software. B. Sterr performed the simulations, analyzed and
visualized the data, and drafted the manuscript. The original manuscript draft
was extensively reviewed and edited by B. Sterr, A. Hrymak, M. Schneider, and
T. Böhlke. Resources were provided by M. Schneider and T. Böhlke. The research
project was supervised by A. Hrymak, M. Schneider and T. Böhlke.

## Appendix

###  A Coated layered materials with singular core

####  A.1 Singularity condition for coated layered materials with singular
core

In the particular cases that the core material $m_{1}$ of a CLM is rigid or
represents a void, the dual or the primal material properties
${\mathbb{K}}_{1}$ or ${\mathbb{M}}_{1}$ approach zero, respectively. Thus,
more generally, we let the material tensor ${\mathbb{A}}_{1}$ approach zero,
and rewrite equation (2.42) as

$\lim_{{\mathbb{A}}_{1}\rightarrow\bm{0}}\bar{{\mathbb{A}}}={\mathbb{A}}_{2}\left({\mathbb{I}}+(1-f_{2})\left(-{\mathbb{I}}+f_{2}\sum_{r=1}^{R}c_{r}{\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{A}}_{2}}(\mbox{\boldmath$n$}_{r}){\mathbb{A}}_{2}\right)^{-1}\right).$
(5.1)

In the following, we study under which conditions the effective material
properties $\bar{{\mathbb{A}}}$ of a rank-$R$ CLM are singular. The effective
material properties $\bar{{\mathbb{A}}}$ are singular precisely if there
exists an effective tensor
$\bar{\mbox{\boldmath$g$}}\in\\{\sf{Sym}(3),{\sf{Sym}_{0}(3)}\\}$, such that

$\bar{{\mathbb{A}}}[\bar{\mbox{\boldmath$g$}}]=\mbox{\boldmath$0$}.$ (5.2)

Using the definition (5.1), equation (5.2) may be equivalently rewritten in
the form

$\left(\sum_{r=1}^{R}c_{r}{\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{A}}_{2}}(\mbox{\boldmath$n$}_{r}){\mathbb{A}}_{2}\right)[\bar{\mbox{\boldmath$g$}}]=\bar{\mbox{\boldmath$g$}}.$
(5.3)

Equation (5.3) is satisfied precisely if $\bar{\mbox{\boldmath$g$}}$ lies in
the intersection of all subspaces ${\mathcal{E}}^{\sf{A}}_{r}$, such that

$\bar{\mbox{\boldmath$g$}}\in\bigcap_{r=1}^{R}{\mathcal{E}}^{\sf{A}}_{r},$
(5.4)

or, equivalently,

${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{I}}}(\mbox{\boldmath$n$}_{r})[\bar{\mbox{\boldmath$g$}}]=\bar{\mbox{\boldmath$g$}},\quad\forall
r\in\\{1,...,R\\}.$ (5.5)

We prove this as follows. First, if $\bar{\mbox{\boldmath$g$}}$ lies in the
section of all subspaces such that (5.4) is satisfied, it follows from
equation (5.5) and the definition of the operator
${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{A}}_{2}}(\mbox{\boldmath$n$}_{r})$
(2.41) that

$\left(\sum_{r=1}^{R}c_{r}{\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{A}}_{2}}(\mbox{\boldmath$n$}_{r}){\mathbb{A}}_{2}\right)[\bar{\mbox{\boldmath$g$}}]=\left(\sum_{r=1}^{R}c_{r}{\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{A}}_{2}}(\mbox{\boldmath$n$}_{r}){\mathbb{A}}_{2}{\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{I}}}(\mbox{\boldmath$n$}_{r})[\bar{\mbox{\boldmath$g$}}]\right)=\sum_{r=1}^{R}c_{r}\bar{\mbox{\boldmath$g$}}=\bar{\mbox{\boldmath$g$}},$
(5.6)

and hence equation (5.3) is satisfied. Second we wish to prove that the
converse is also true, i.e., if equation (5.3) is satisfied, equation (5.4) is
satisfied as well. Let ${\mathbb{T}}_{{\mathbb{A}}_{2}}$ denote a convex
combination of the projectors
${\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{A}}_{2}}(\mbox{\boldmath$n$}_{r}){\mathbb{A}}_{2}$,
i.e.,

${\mathbb{T}}_{{\mathbb{A}}_{2}}=\left(\sum_{r=1}^{R}c_{r}{\mathchoice{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\displaystyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\textstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptstyle\mathbb{L}$}}}{\raisebox{0.0pt}{\scalebox{1.0}[-1.0]{$\scriptscriptstyle\mathbb{L}$}}}}_{{\mathbb{A}}_{2}}(\mbox{\boldmath$n$}_{r}){\mathbb{A}}_{2}\right),$
(5.7)

with

$c_{r}>0\quad\text{and}\quad\sum_{r=1}^{R}c_{r}=1.$ (5.8)

Then, fixed points of the operator ${\mathbb{T}}_{{\mathbb{A}}_{2}}$ lie in
the set [94, Lemma 1.4]

$\mathcal{E}=\bigcap_{r=1}^{R}{\mathcal{E}}^{\sf{A}}_{r}.$ (5.9)

Therefore, any eigenvector of the operator ${\mathbb{T}}_{{\mathbb{A}}_{2}}$
corresponding to the eigenvalue one must lie in the intersection of the
subspaces ${\mathcal{E}}^{\sf{A}}_{r}$, which we denote with $\mathcal{E}$. In
particular, if equation (5.3) is satisfied, the element
$\bar{\mbox{\boldmath$g$}}$ lies in $\mathcal{E}$, i.e.,

${\mathbb{T}}_{{\mathbb{A}}_{2}}[\bar{\mbox{\boldmath$g$}}]=\bar{\mbox{\boldmath$g$}},\quad\implies\quad\bar{\mbox{\boldmath$g$}}\in\mathcal{E}$
(5.10)

and equation (5.4) is satisfied as well.

With the result (5.10) at hand, we observe that a necessary and sufficient
condition for the regularity of a rank-R CLM is that the space $\mathcal{E}$
is trivial, i.e.,

$\mathcal{E}=\\{\mbox{\boldmath$0$}\\}.$ (5.11)

We refer to equation (5.11) as non-singularity condition in the following. In
the context of linear elastic and linearly viscous materials, we derived the
non-singularity condition (5.11) for the spaces ${\mathcal{E}}^{\sf{A}}_{r}$,
such that

${\mathcal{E}}^{\sf{A}}_{r}\subset\sf{Sym}(3),\quad\forall