# A reduced order model for geometrically parameterized two-scale simulations
of elasto-plastic microstructures under large deformations

Theron Guo<EMAIL_ADDRESS>Ondřej Rokoš<EMAIL_ADDRESS>Karen Veroy
<EMAIL_ADDRESS>Centre for Analysis, Scientific Computing and Applications,
Eindhoven University of Technology, 5612 AZ Eindhoven, The Netherlands
Mechanics of Materials, Eindhoven University of Technology, 5612 AZ Eindhoven,
The Netherlands  Institute for Complex Molecular Systems, Eindhoven University
of Technology, 5612 AZ Eindhoven, The Netherlands

###### Abstract

In recent years, there has been a growing interest in understanding complex
microstructures and their effect on macroscopic properties. In general, it is
difficult to derive an effective constitutive law for such microstructures
with reasonable accuracy and meaningful parameters. One numerical approach to
bridge the scales is computational homogenization, in which a microscopic
problem is solved at every macroscopic point, essentially replacing the
effective constitutive model. Such approaches are, however, computationally
expensive and typically infeasible in multi-query contexts such as
optimization and material design. To render these analyses tractable,
surrogate models that can accurately approximate and accelerate the
microscopic problem over a large design space of shapes, material and loading
parameters are required. In previous works, such models were constructed in a
data-driven manner using methods such as Neural Networks (NN) or Gaussian
Process Regression (GPR). However, these approaches currently suffer from
issues, such as need for large amounts of training data, lack of physics, and
considerable extrapolation errors. In this work, we develop a reduced order
model based on Proper Orthogonal Decomposition (POD), Empirical Cubature
Method (ECM) and a geometrical transformation method with the following key
features: (i) large shape variations of the microstructure are captured, (ii)
only relatively small amounts of training data are necessary, and (iii) highly
non-linear history-dependent behaviors are treated. The proposed framework is
tested and examined in two numerical examples, involving two scales and large
geometrical variations. In both cases, high speed-ups and accuracies are
achieved while observing good extrapolation behavior.

###### keywords:

Reduced order modelling , proper orthogonal decomposition , computational
homogenization , hyperreduction , empirical cubature method , geometrical
transformation

capbtabboxtable[][] [table]capposition=top

##  1 Introduction

Driven by advances in additive manufacturing and tailorable effective
properties of metamaterials, there has been a growing interest in
understanding structure-property relationships of complex microstructures.
These microstructures can typically be described by a few shape parameters,
leading to distinct types of effective behavior. To investigate such
structure-property relations and to find the optimal shape for a given
application, simulations are often considered. These simulations are in
general computationally expensive or even intractable for direct numerical
simulation, especially for large engineering applications, since considerably
fine meshes are required to capture the complex microstructural geometry. By
employing multi-scale methods [1, 2] or domain decomposition methods [3, 4],
such large-scale problems can be separated into many smaller subproblems, thus
rendering them amenable for efficient numerical simulation.

If scale separation is assumed, i.e., when the length scale of the typical
microstructural features is much smaller than that of the macrostructure,
first-order computational homogenization can be employed. Here, the behavior
of the microstructure dictates the (average) constitutive behavior of an
effective macrostructural continuum model. By defining a Representative Volume
Element (RVE) which models the fine-scale geometry of the microstructure in
full detail, a coarse-grained representation of the macrostructure with a much
coarser discretization can be assumed at the macroscale. At every macroscopic
integration point, the macroscopic strain is used to specify a microscopic
boundary value problem which, after solution, returns the effective stress and
stiffness. Since a partial differential equation (PDE) needs to be solved at
every macroscopic Gauss integration point, this methodology is still
computationally expensive, and efficient ways for its solution are needed.

Several approaches to tackle this problem have been reported in the
literature. For instance, the Fast Fourier Transform (FFT) [5, 6] allows to
directly use pixelized images of real microstructures, thus avoiding costly
meshing of complex microstructures. The solution can further be accelerated by
the (nonuniform) Transformation Field Analysis (see, e.g., [7, 8]), or Self-
consistent Clustering Analysis [9, 10]. One disadvantage of FFT is that
geometrical parameterizations of the RVE cannot be directly treated and,
hence, sensitivities for material optimization cannot be directly computed.
Another class of methods aims at solving the microscopic problem via the
Finite Element (FE) method, resulting in a multi-scale formulation referred to
as FE2 [11, 12, 13]. By directly solving the microscopic PDE, material or
shape parameterizations can be considered in a straightforward manner, making
the approach more suitable for inverse problems and optimization. To speed up
the microstructural simulation, proper orthogonal decomposition (POD) [14, 15]
can be utilized to find a reduced set of basis functions; the method then
computes the Galerkin projection of the solution onto the space spanned by the
snapshots. Although POD generally requires many full-order solves for
training, it typically works well for all input parameters. In the context of
first-order homogenization, POD was first applied in Yvonnet et al. [16] for a
hyper-elastic RVE, and later explored in Radermacher et al. [17] for an
elasto-plastic RVE under small strains. However, due to the non-linearities of
the microscopic problem, the speed-ups were limited since the global force
vector and stiffness matrix must be assembled by full integration in every
microscopic Newton iteration. To address this issue, a further reduction
called hyperreduction is required, which aims at finding an efficient way of
assembling microstructural force and stiffness quantities. Notable
hyperreduction methods are empirical interpolation method (EIM) (see, e.g.,
[18]), a variant of EIM called discrete empirical interpolation method (DEIM)
(see, e.g., [19]), energy-based mesh sampling and weighting [20], reduced
integration domain [21], empirical quadrature procedure [22], and empirical
cubature method (ECM) [23]. EIM and DEIM interpolate the non-linear integrand
of the global force vector such that the integrals can be pre-computed. In
[24, 25], DEIM was used successfully to accelerate the solution of the
microscopic PDE. However, these works only discussed the solution of the
microscopic PDE and did not derive the effective stress and stiffness
quantities required for the macroscopic problem. A possible disadvantage of
EIM and DEIM is that they lead to non-symmetric tangent matrices, which might
result in convergence issues, observed in, for instance, [25, 26]. The
remainder of the above-mentioned hyperreduction methods aim at approximating
the integrals by finding a subset of integration points with corresponding
weights among the set of all integration points used in the formulation of the
microstructural PDE. This has the advantage that the stiffness matrix is
always symmetric, ensuring a good convergence of the microscopic problem.
Example applications of hyperreduction methods have been successfully employed
in two-scale simulations in [27], where an elasto-plastic composite RVE under
large deformations was considered. In [28], a damage model for a composite RVE
under small deformations was shown. While both works obtained accurate results
and successfully accelerated the forward simulations of a two-scale problem,
such formulations were limited to fixed microstructures only, i.e., did not
account for possible parameterizations. In order to allow for optimization of
microstructures, the surrogate model needs to be extended to a wide range of
different design parameters (including geometrical as well as material).

This work aims to address this gap, by developing a hyper-reduced surrogate
model for geometrically parameterized microstructures, to enable (shape)
sensitivity analysis and optimization of materials. Furthermore, we intend to
provide a detailed analysis of the reduced RVE problem for arbitrary loading
paths and geometries and elucidate possible issues due to reduction. Our main
contributions are:

  1. 1.

development of a hyper-reduced POD model for a family of geometrically
parameterized microstructures, by employing a geometrical transformation
method [29] and by extending the ECM algorithm to geometrical parameters,

  2. 2.

consistent derivation of the effective stress and stiffness of the hyper-
reduced model,

  3. 3.

an empirical analysis of the accuracy of the surrogate model for elasto-
plastic RVEs under large deformations for different geometries and loading
conditions,

  4. 4.

a quantitative comparison of a two-scale example with continuous change in
microstructural heterogeneities.

The remainder of this paper is organized as follows. In Section 2, the
microscopic problem arising in first-order computational homogenization is
briefly summarized, together with the computation of the effective quantities.
Section 3 covers in-depth development of the reduced order model with
particular focus on the empirical cubature method for geometrically
parameterized microstructures, and includes a detailed derivation of the
effective stress and stiffness. In Section 4, the proposed method is examined
and tested in detail, first for a single RVE and then also for a full two-
scale problem. Finally, a summary on the findings and concluding remarks are
given in Section 5.

In this work, the following notational convention is adopted. Italic bold
symbols are used for coordinates $\bm{X}$ and vectorial or tensorial fields,
such as the displacement field $\bm{u}$ or stress field $\bm{P}$. Upright bold
symbols are used for algebraic vectors and matrices, such as the global
stiffness matrix $\mathbf{K}$ or the coefficients of the discretized
displacement field $\mathbf{u}$. A field quantity $\bm{u}$ for given
parameters $\bm{\mu}$ is denoted as $\bm{u}(\bm{X};\bm{\mu})$. Given second-
order tensors $\bm{A}$ and $\bm{B}$, fourth-order tensor $\bm{C}$, and vector
$\bm{v}$, the following operations are used: $\bm{A}\bm{B}=A_{ij}B_{jk}$,
$\bm{A}:\bm{B}=A_{ij}B_{ij}$, $\bm{A}:\bm{C}:\bm{B}=A_{ij}C_{ijkl}B_{kl}$ and
$\bm{A}\bm{v}=A_{ij}v_{j}$, where the Einstein summation convention is
implied.

##  2 Formulation Of The Microscopic Problem

In multiscale schemes based on first-order homogenization, the macroscopic
constitutive model is replaced by a microscopic PDE which is defined on an
RVE. By prescribing the macroscopic deformation gradient on the microscale,
the PDE can be solved and an effective stress and stiffness is returned to the
macroscopic solver, see Fig. 1. For applications such as microstructure
optimization, it is reasonable to additionally introduce a parameterization of
the RVE in order to compute the sensitivities with respect to design
variables. The microscopic boundary value problem is formulated below on a
parameterized domain, as is usually the case in shape optimization. For
brevity, the dependence on the macroscopic coordinates is omitted and a fixed
macroscopic material point is assumed unless otherwise specified.

###  2.1 Boundary Value Problem

Figure 1: Two-scale problem based on first-order homogenization. At every
macroscopic point, a microscopic simulation is defined through deformation
gradient $\bar{\bm{F}}$ and shape parameters $\bm{\mu}$, and solved to obtain
an effective stress $\bar{\bm{P}}$ and stiffness $\bar{\bm{A}}$. For different
macroscopic points, different parameterized microstructures can be considered
through $\bm{\mu}$. As an example of a family of geometrically parameterized
microstructures, a parent domain with a circular inclusion $\Omega^{\rm{p}}$
(center), can be mapped onto parameterized domains $\Omega^{\bm{\mu}_{1}}$
(left) and $\Omega^{\bm{\mu}_{2}}$ (right) with mapping
$\bm{\Phi}_{\bm{\mu}_{1}}$ and $\bm{\Phi}_{\bm{\mu}_{2}}$.

Consider a family of domains $\Omega^{\bm{\mu}}\subset\mathbb{R}^{d}$ with
space dimension $d=2,3$ the space dimension, parameterized by geometrical
parameters $\bm{\mu}$, and spanned by a position vector
$\bm{X}^{\bm{\mu}}\in\Omega^{\bm{\mu}}$. In Fig. 1, an example parent domain
with a circular inclusion $\Omega^{\rm{p}}$ is geometrically parameterized and
mapped to two distinct parameterized domains with elliptical inclusions,
$\Omega^{\bm{\mu}_{1}}$ and $\Omega^{\bm{\mu}_{2}}$. The volume and the
topology of the domain $|\Omega^{\bm{\mu}}|$ are assumed to remain fixed for
all parameters (the outer boundaries of the RVE domain are fixed while the
shape of the interior geometry can change). With the assumption of scale
separation between macro- and microscale, the microscopic displacement field
on the parameterized domain $\bm{u}(\bm{X}^{\bm{\mu}})$ can be written as the
summation of a mean field $\bar{\bm{u}}(\bm{X}^{\bm{\mu}})$ and a fluctuation
field $\bm{w}(\bm{X}^{\bm{\mu}})$, i.e.,
$\bm{u}(\bm{X}^{\bm{\mu}})=\bar{\bm{u}}(\bm{X}^{\bm{\mu}})+\bm{w}(\bm{X}^{\bm{\mu}})$.
The mean field is fully specified through
$\bar{\bm{u}}(\bm{X}^{\bm{\mu}})\coloneqq(\bar{\bm{F}}-\bm{I})\bm{X}^{\bm{\mu}}$,
where $\bar{\bm{F}}$ is the macroscopic deformation gradient tensor and
$\bm{I}$ is the identity tensor. The total deformation gradient tensor
$\bm{F}$ is defined as

$\displaystyle\bm{F}(\bm{w}(\bm{X}^{\bm{\mu}}))\coloneqq\dfrac{\partial\bm{u}}{\partial\bm{X}^{\bm{\mu}}}=\bar{\bm{F}}+\dfrac{\partial\bm{w}}{\partial\bm{X}^{\bm{\mu}}}.$
(1)

The governing microscopic PDE is given as

$\displaystyle\begin{aligned}
\operatorname{Div}\bm{P}^{T}(\bm{F}(\bm{w}(\bm{X}^{\bm{\mu}})))&=\bm{0}&&\text{
on }\Omega^{\bm{\mu}},\\\ \bm{w}&\text{ periodic }\quad&&\text{ on
}\partial\Omega^{\bm{\mu}},\end{aligned}$ (2)

where $\operatorname{Div}(\bullet)$ is the divergence operator with respect to
$\bm{X}^{\bm{\mu}}$ and $\bm{P}$ denotes the second-order first Piola-
Kirchhoff (1PK) stress tensor. No constitutive model is specified at this
point, although we assume that the stress $\bm{P}$ is a non-linear function of
the deformation gradient $\bm{F}$ (or its history). The weak form of the
problem is then: given the macroscopic deformation gradient $\bar{\bm{F}}$,
find the fluctuation field
$\bm{w}^{*}\in\mathcal{V}\coloneqq\\{\bm{v}\in(H^{1}(\Omega^{\bm{\mu}}))^{d}\
|\ \bm{v}\text{ periodic on }\partial\Omega^{\bm{\mu}}\\}$ that fulfills

$\displaystyle
G(\bm{w})\coloneqq\int_{\Omega^{\bm{\mu}}}\frac{\partial\delta\bm{w}}{\partial\bm{X}^{\bm{\mu}}}:\bm{P}\left(\bar{\bm{F}}+\frac{\partial\bm{w}}{\partial\bm{X}^{\bm{\mu}}}\right)d\bm{X}^{\bm{\mu}}\stackrel{{\scriptstyle!}}{{=}}0,\qquad\forall\delta\bm{w}\in\mathcal{V},$
(3)

where the integral bounds depend on the parameters $\bm{\mu}$, $\delta\bm{w}$
denotes a test function, and $H^{1}(\Omega^{\bm{\mu}})$ is a Hilbert space
with square integrable functions and square integrable derivatives. The inner
product in $\mathcal{V}$ is defined as

$\displaystyle(\bm{u},\bm{v})_{\mathcal{V}}\coloneqq\int_{\Omega^{\bm{\mu}}}\left(\bm{u}\cdot\bm{v}+\frac{\partial\bm{u}}{\partial\bm{X}^{\bm{\mu}}}:\frac{\partial\bm{v}}{\partial\bm{X}^{\bm{\mu}}}\right)d\bm{X}^{\bm{\mu}}.$
(4)

From Eq. 3, it is apparent that the macroscopic deformation gradient
$\bar{\bm{F}}$ represents the external loading, while the fluctuation
displacement field $\bm{w}$ balances the system. To simplify the problem in
Eq. 3 and remove the parameter dependence of the integral bounds, a parent
domain $\Omega^{\rm{p}}$ is defined. To this end, we assume that there exists
a parameter-dependent diffeomorphism
$\bm{\Phi}_{\bm{\mu}}:\Omega^{\rm{p}}\rightarrow\Omega^{\bm{\mu}},\bm{X}^{\rm{p}}\mapsto\bm{X}^{\bm{\mu}}$,
see Fig. 1. Using integration by substitution, the problem of Eq. 3 can be
restated as follows: given the macroscopic deformation gradient
$\bar{\bm{F}}$, find
$\bm{w}^{*\rm{p}}\in\mathcal{V}^{\rm{p}}\coloneqq\\{\bm{v}\in(H^{1}(\Omega^{\rm{p}}))^{d}\
|\ \bm{v}\text{ periodic on }\partial\Omega^{\rm{p}}\\}$ that fulfills

$\displaystyle
G^{\rm{p}}(\bm{w}^{\rm{p}})\coloneqq\int_{\Omega^{\rm{p}}}\left(\frac{\partial\delta\bm{w}^{\rm{p}}}{\partial\bm{X}^{\rm{p}}}\bm{F}_{\bm{\mu}}^{-1}\right):\bm{P}\left(\bar{\bm{F}}+\frac{\partial\bm{w}^{\rm{p}}}{\partial\bm{X}^{\rm{p}}}\bm{F}_{\bm{\mu}}^{-1}\right)\left|\det\bm{F}_{\bm{\mu}}\right|d\bm{X}^{\rm{p}}=0,\qquad\forall\delta\bm{w}^{\rm{p}}\in\mathcal{V}^{\rm{p}},$
(5)

with the transformation gradient
$\bm{F}_{\bm{\mu}}\coloneqq\dfrac{\partial\bm{\Phi_{\bm{\mu}}}}{\partial\bm{X}^{\rm{p}}}$
and $d\bm{X}^{\bm{\mu}}=\left|\det\bm{F}_{\bm{\mu}}\right|d\bm{X}^{\rm{p}}$.
The superscript p is used to denote quantities pertinent to the parent domain,
e.g.,
$\bm{w}(\bm{X}^{\bm{\mu}})=(\bm{w}\circ\bm{\Phi}_{\bm{\mu}})(\bm{X}^{\rm{p}})=\bm{w}^{\rm{p}}(\bm{X}^{\rm{p}})$.
To iteratively solve the non-linear problem in Eq. 5, a linearization using
the Gateaux derivative around the current state $\bm{w}^{\rm{p}}$ in direction
$\Delta\bm{w}^{\rm{p}}\in\mathcal{V}^{\rm{p}}$ is required and can be written
as,

$\displaystyle\begin{aligned} \left.\frac{\partial
G^{\rm{p}}(\bm{w}^{\rm{p}}+\tau\Delta\bm{w}^{\rm{p}})}{\partial\tau}\right|_{\tau=0}&=\int_{\Omega^{\rm{p}}}\left(\frac{\partial\delta\bm{w}^{\rm{p}}}{\partial\bm{X}^{\rm{p}}}\bm{F}_{\bm{\mu}}^{-1}\right):\bm{A}\left(\bar{\bm{F}}+\frac{\partial\bm{w}^{\rm{p}}}{\partial\bm{X}^{\rm{p}}}\bm{F}_{\bm{\mu}}^{-1}\right):\left(\frac{\partial\Delta\bm{w}^{\rm{p}}}{\partial\bm{X}^{\rm{p}}}\bm{F}_{\bm{\mu}}^{-1}\right)\left|\det\bm{F}_{\bm{\mu}}\right|d\bm{X}^{\rm{p}},\end{aligned}$
(6)

where $\bm{A}\coloneqq\dfrac{\partial\bm{P}}{\partial\bm{F}}$ is the fourth-
order stiffness tensor. Once the transformation map $\bm{\Phi}_{\bm{\mu}}$ is
known, Eq. 3 can be solved on the parent domain using Eqs. 5 and 6. Further
details on how to find these transformations for a range of geometrical
parameters are provided in Section 3.4.

By employing a finite element discretization for
$\bm{w}^{\rm{p}}\approx\bm{w}^{\rm{p}}_{h}\in\mathcal{V}_{h}^{\rm{p}}\subset\mathcal{V}^{\rm{p}}$,
with $\dim\mathcal{V}_{h}^{\rm{p}}=\mathcal{N}$, the number of degrees of
freedom of the dicretization, the internal force vector
$\mathbf{f}\in\mathbb{R}^{\mathcal{N}}$ and global stiffness matrix
$\mathbf{K}\in\mathbb{R}^{\mathcal{N}\times\mathcal{N}}$ can be derived from
Eqs. 5 and 6, resulting in the following non-linear system of equations

$\displaystyle\mathbf{f}(\mathbf{w})=\bm{0},$ (7)

where $\mathbf{w}\in\mathbb{R}^{\mathcal{N}}$ is the column vector of unknown
coefficients of the discretized fluctuation field. This problem can be solved
with the Newton method, i.e.,

$\displaystyle\begin{aligned}
\mathbf{K}(\mathbf{w}^{m})\Delta\mathbf{w}&=-\mathbf{f}(\mathbf{w}^{m}),\\\
\mathbf{w}^{m+1}&=\mathbf{w}^{m}+\Delta\mathbf{w},\end{aligned}$ (8)

where $m$ is the number of Newton iteration and Eq. 8 is repeated until
$||\mathbf{f}(\mathbf{w}^{m})||_{2}\leq\varepsilon_{\text{newton}}$ with
$\varepsilon_{\text{newton}}$ a user-defined tolerance.

####  2.1.1 Effective Quantities

For conciseness of notation, the following abbreviations are introduced to
denote quantities after the solution $\bm{w}^{*\rm{p}}$ has been obtained:

$\displaystyle\bm{P}^{*\rm{p}}$
$\displaystyle\coloneqq\bm{P}\left(\bar{\bm{F}}+\frac{\partial\bm{w}^{*\rm{p}}}{\partial\bm{X}^{\rm{p}}}\bm{F}_{\bm{\mu}}^{-1}\right),$
(9) $\displaystyle\bm{A}^{*\rm{p}}$
$\displaystyle\coloneqq\bm{A}\left(\bar{\bm{F}}+\frac{\partial\bm{w}^{*\rm{p}}}{\partial\bm{X}^{\rm{p}}}\bm{F}_{\bm{\mu}}^{-1}\right).$
(10)

Upon obtaining solution $\bm{w}^{*\rm{p}}$ from Eq. 5, the effective stress is
computed as

$\displaystyle\bar{\bm{P}}$
$\displaystyle\coloneqq|\Omega^{\rm{p}}|^{-1}\int_{\Omega^{\rm{p}}}\bm{P}^{*\rm{p}}|\det\bm{F}_{\bm{\mu}}|d\bm{X}^{\rm{p}},$
(11)

and the effective stiffness (in index notation) as

$\displaystyle\begin{aligned}
\bar{A}_{ijkl}\coloneqq&\frac{\partial\bar{P}_{ij}}{\partial\bar{F}_{kl}}\\\
=&|\Omega^{\rm{p}}|^{-1}\frac{\partial}{\partial\bar{F}_{kl}}\int_{\Omega^{\rm{p}}}P^{*\rm{p}}_{ij}|\det\bm{F}_{\bm{\mu}}|d\bm{X}^{\rm{p}}\\\
=&|\Omega^{\rm{p}}|^{-1}\int_{\Omega^{\rm{p}}}A^{*\rm{p}}_{ijmn}\left(\mathbb{I}_{mnkl}+\frac{\partial}{\partial\bar{F}_{kl}}\left(\frac{\partial
w_{m}^{*}}{\partial
X_{r}}\right)\left(F_{\bm{\mu}}^{-1}\right)_{rn}\right)|\det\bm{F}_{\bm{\mu}}|d\bm{X}^{\rm{p}},\end{aligned}$
(12)

where $\mathbb{I}_{mnkl}\coloneqq\delta_{mk}\delta_{nl}$ is the fourth-order
identity tensor. To determine
$\dfrac{\partial}{\partial\bar{F}_{kl}}\left(\dfrac{\partial
w_{m}^{*}}{\partial X_{r}}\right)$, Eq. 5 is differentiated with respect to
$\bar{\bm{F}}$. For one particular component $\bar{F}_{kl}$ (where the indices
$k$ and $l$ are assumed to be temporarily fixed), the differentiation yields

$\displaystyle\int_{\Omega^{\rm{p}}}\left(\frac{\partial\delta\bm{w}^{\rm{p}}}{\partial\bm{X}^{\rm{p}}}\bm{F}_{\bm{\mu}}^{-1}\right):\bm{A}^{*\rm{p}}:\left(\frac{\partial\bm{q}_{kl}}{\partial\bm{X}^{\rm{p}}}\bm{F}_{\bm{\mu}}^{-1}\right)\left|\det\bm{F}_{\bm{\mu}}\right|d\bm{X}^{\rm{p}}=-\left(\int_{\Omega^{\rm{p}}}\left(\frac{\partial\delta\bm{w}^{\rm{p}}}{\partial\bm{X}^{\rm{p}}}\bm{F}_{\bm{\mu}}^{-1}\right):\bm{A}^{*\rm{p}}\left|\det\bm{F}_{\bm{\mu}}\right|d\bm{X}^{\rm{p}}\right):\mathbb{E}_{kl},$
(13)

where a new auxiliary vector field
$\bm{q}_{kl}\coloneqq\dfrac{\partial\bm{w}^{*\rm{p}}}{\partial\bar{F}_{kl}}\in\mathcal{V}^{\rm{p}}$
has been defined (reflecting the sensitivity of the microfluctuation field
with respect to the change of the applied macroscopic loading), and
$\mathbb{E}_{kl}\in\mathbb{R}^{d\times d}$ is a second order tensor with all
entries zero, except for the $kl$-th entry which is 1. The linear tangent
problem of Eq. 13 is then solved for all combinations $k,l=1,...,d$ to obtain
$\bm{q}_{kl}$ for each component of $\bar{\bm{F}}$.

For shape optimization, the sensitivities of the effective stress
$\bar{\bm{P}}$ with respect to the geometrical parameters $\bm{\mu}$ are
required and are computed as follows (in index notation)

$\displaystyle\frac{\partial\bar{P}_{ij}}{\partial\mu_{k}}=|\Omega^{\rm{p}}|^{-1}\int_{\Omega^{\rm{p}}}\left(A^{*\rm{p}}_{ijmn}\left(\frac{\partial}{\partial\mu_{k}}\left(\frac{\partial
w_{m}^{*}}{\partial
X_{r}}\right)\left(F_{\bm{\mu}}^{-1}\right)_{rn}+\frac{\partial
w_{m}^{*}}{\partial
X_{r}}\frac{\partial\left(F_{\bm{\mu}}^{-1}\right)_{rn}}{\partial\mu_{k}}\right)|\det\bm{F}_{\bm{\mu}}|+P^{*\rm{p}}_{ij}\frac{\partial|\det\bm{F}_{\bm{\mu}}|}{\partial\mu_{k}}\right)d\bm{X}^{\rm{p}}.$
(14)

The left-hand expression is complicated due to the derivatives of
$\bm{F}_{\bm{\mu}}^{-1}$ and $|\det\bm{F}_{\bm{\mu}}|$. If the effective
stress $\bar{\bm{P}}$ can be assumed to vary smoothly with the parameters
$\bm{\mu}$, which may be a reasonable assumption for smoothly varying shapes
using, for instance, splines, finite differences can be used to approximate
these sensitivities.

##  3 Surrogate Modelling

Since the microscopic problem has to be solved at every macroscopic quadrature
point, the solution of the microscopic PDE must be efficient. Solving it
directly using FE is in general too computationally expensive, and, hence, the
microscopic solver must be accelerated. In this section, a surrogate model for
the geometrically parameterized microscopic PDE is developed by employing the
reduced basis method (RBM) [14] to reduce the number of degrees of freedom and
the empirical cubature method (ECM) [23] to reduce the number of quadrature
points. The key idea is to construct the surrogate model on the parent domain
$\Omega^{\rm{p}}$, adapt it to each geometry $\Omega^{\bm{\mu}}$, and then
solve the reduced problem.

###  3.1 Reduced Basis Method

For complex problems and geometries, typically a fine mesh is required for FE,
leading to a high-dimensional solution space $\mathcal{V}_{h}^{\rm{p}}$ for
the fluctuation displacement field $\bm{w}^{\rm{p}}$ with
$\dim\mathcal{V}_{h}^{\rm{p}}=\mathcal{N}$. The idea of the RBM is to
approximate the field with global parameter-independent basis functions and
parameter-dependent coefficients, i.e.,

$\displaystyle\bm{w}^{\rm{p}}(\bm{X}^{\rm{p}};\bar{\bm{F}},\bm{\mu})\approx\sum_{n=1}^{N}a_{n}(\bar{\bm{F}},\bm{\mu})\bm{\phi}_{n}(\bm{X}^{\rm{p}}),$
(15)

where $N$ is the number of basis functions, ideally much smaller than the
dimension of the FE space, i.e., $N\ll\mathcal{N}$. The basis functions,
$\\{\bm{\phi}_{n}\\}_{n=1}^{N}$, span a subset of $\mathcal{V}^{\rm{p}}_{h}$
and can be obtained by applying proper orthogonal decomposition (POD) on a set
of pre-computed full solutions for different parameter values. Additionally,
they are orthonormal with respect to $\mathcal{V}^{\rm{p}}$, i.e.,

$\displaystyle(\bm{\phi}_{m},\bm{\phi}_{n})_{\mathcal{V}^{\rm{p}}}=\delta_{mn},$
(16)

where $\delta_{mn}$ denotes the Kronecker delta. By utilizing the POD space
for both the trial and test space and inserting $\bm{w}^{\rm{p}}$ from Eq. 15
into Eqs. 5 and 6, the components for the reduced internal force vector
$\mathbf{f}^{\text{POD}}\in\mathbb{R}^{N}$ and reduced global stiffness matrix
$\mathbf{K}^{\text{POD}}\in\mathbb{R}^{N\times N}$ can be derived as

$\displaystyle f^{\text{POD}}_{i}(\mathbf{a})$
$\displaystyle\coloneqq\int_{\Omega^{\rm{p}}}\left(\frac{\partial\bm{\phi}_{i}}{\partial\bm{X}^{\rm{p}}}\bm{F}_{\bm{\mu}}^{-1}\right):\bm{P}\left(\bar{\bm{F}}+\left(\sum_{n=1}^{N}a_{n}\frac{\partial\bm{\phi}_{n}}{\partial\bm{X}^{\rm{p}}}\right)\bm{F}_{\bm{\mu}}^{-1}\right)\left|\det\bm{F}_{\bm{\mu}}\right|d\bm{X}^{\rm{p}},$
(17) $\displaystyle K^{\text{POD}}_{ij}(\mathbf{a})$
$\displaystyle\coloneqq\int_{\Omega^{\rm{p}}}\left(\frac{\partial\bm{\phi}_{i}}{\partial\bm{X}^{\rm{p}}}\bm{F}_{\bm{\mu}}^{-1}\right):\bm{A}\left(\bar{\bm{F}}+\left(\sum_{n=1}^{N}a_{n}\frac{\partial\bm{\phi}_{n}}{\partial\bm{X}^{\rm{p}}}\right)\bm{F}_{\bm{\mu}}^{-1}\right):\left(\frac{\partial\bm{\phi}_{j}}{\partial\bm{X}^{\rm{p}}}\bm{F}_{\bm{\mu}}^{-1}\right)\left|\det\bm{F}_{\bm{\mu}}\right|d\bm{X}^{\rm{p}},$
(18)

where $\mathbf{a}=[a_{1},\dots,a_{N}]^{T}$ is the column vector of unknown
coefficients to be solved for, and $i,j=1,\dots,N$ span over all basis
functions. Analogously to Eqs. 7 and 8, the resulting non-linear system of
equations

$\displaystyle\mathbf{f}^{\text{POD}}(\mathbf{a})=\bm{0}$ (19)

can be solved using Newton method:

$\displaystyle\begin{aligned}
\mathbf{K}^{\text{POD}}(\mathbf{a}^{m})\Delta\mathbf{a}&=-\mathbf{f}^{\text{POD}}(\mathbf{a}^{m}),\\\
\mathbf{a}^{m+1}&=\mathbf{a}^{m}+\Delta\mathbf{a}.\end{aligned}$ (20)

###  3.2 Empirical Cubature Method

Even though the solution field and linear system of equations have been
reduced to dimension $N\ll\mathcal{N}$, computing the components of the force
vector in Eq. 17 and global stiffness matrix in Eq. 18 still requires
integrating over the RVE. For the full integration, a numerical quadrature
rule (usually based on Gauss quadrature) with integration points and
corresponding weights $\\{(\hat{\bm{X}}_{q},\hat{w}_{q})\\}_{q=1}^{\hat{Q}}$,
where $\hat{Q}$ is the total number of integration points, is employed, i.e.,

$\displaystyle
f^{\text{POD}}_{i}(\mathbf{a})\approx\sum_{q=1}^{\hat{Q}}\hat{w}_{q}\left.\left[\left(\frac{\partial\bm{\phi}_{i}}{\partial\bm{X}^{\rm{p}}}\bm{F}_{\bm{\mu}}^{-1}\right):\bm{P}\left(\bar{\bm{F}}+\left(\sum_{n=1}^{N}a_{n}\frac{\partial\bm{\phi}_{n}}{\partial\bm{X}^{\rm{p}}}\right)\bm{F}_{\bm{\mu}}^{-1}\right)\left|\det\bm{F}_{\bm{\mu}}\right|\right]\right|_{\hat{\bm{X}}_{q}},$
(21)

for $i=1,\dots,N$. For a fine mesh, $\hat{Q}$ is very large and thus
evaluating Eq. 21 leads to high computational costs. To address this issue, we
employ the Empirical Cubature Method (ECM), which was proposed in Hernández et
al. [23] for a fixed geometry, and extend it to parameterized geometries.

The idea of ECM is to find a subset of points
$\\{\bm{X}_{q}\\}_{q=1}^{Q}\subset\\{\hat{\bm{X}}_{q}\\}_{q=1}^{\hat{Q}}$ with
$Q\ll\hat{Q}$ among the set of all integration points with corresponding
weights $\\{w_{q}\\}_{q=1}^{Q}$ that approximates Eq. 21 up to a user-defined
error $\varepsilon$. To find such a subset that approximates Eq. 21 well for
all admissible geometrical parameters $\bm{\mu}$, Eq. 21 is first rewritten as

$\displaystyle\begin{aligned}
f^{\text{POD}}_{i}(\mathbf{a})&=\sum_{q=1}^{\hat{Q}}\hat{w}_{q}\left.\left[\frac{\partial\bm{\phi}_{i}}{\partial\bm{X}^{\rm{p}}}:\underbrace{\left(\bm{P}\left(\bar{\bm{F}}+\left(\sum_{n=1}^{N}a_{n}\frac{\partial\bm{\phi}_{n}}{\partial\bm{X}^{\rm{p}}}\right)\bm{F}_{\bm{\mu}}^{-1}\right)\bm{F}_{\bm{\mu}}^{-T}\left|\det\bm{F}_{\bm{\mu}}\right|\right)}_{\bm{W}\coloneqq}\right]\right|_{\hat{\bm{X}}_{q}},\\\
&=\sum_{q=1}^{\hat{Q}}\hat{w}_{q}\left.\left[\frac{\partial\bm{\phi}_{i}}{\partial\bm{X}^{\rm{p}}}:\bm{W}(\bm{X}^{\rm{p}};\bar{\bm{F}},\bm{\mu})\right]\right|_{\hat{\bm{X}}_{q}},\end{aligned}$
(22)

where the weighted stress $\bm{W}$ is defined. To remove the parameter
dependence of the integrand in Eq. 22, the weighted stress is approximated by
another reduced basis, i.e.,

$\displaystyle\bm{W}(\bm{X}^{\rm{p}};\bar{\bm{F}},\bm{\mu})\approx\sum_{l=1}^{L}\alpha_{l}(\bar{\bm{F}},\bm{\mu})\bm{B}_{l}(\bm{X}^{\rm{p}}),$
(23)

where $\\{\bm{B}_{l}\\}_{l=1}^{L}$ is a set of $L$ basis functions obtained
using POD, which are orthonormal with respect to $L^{2}(\Omega)$, i.e.,

$\displaystyle\int_{\Omega}\bm{B}_{m}:\bm{B}_{n}d\bm{X}^{\rm{p}}=\delta_{mn}.$
(24)

Inserting Eq. 23 into Eq. 22 and rearranging yields,

$\displaystyle
f^{\text{POD}}_{i}(\mathbf{a})\approx\sum_{l=1}^{L}\alpha_{l}(\bar{\bm{F}},\bm{\mu})\sum_{q=1}^{\hat{Q}}\hat{w}_{q}\left.\left[\frac{\partial\bm{\phi}_{i}}{\partial\bm{X}}:\bm{B}_{l}\right]\right|_{\hat{\bm{X}}_{q}},\qquad
i=1,\dots,N.$ (25)

Since Eq. 25 should be accurate for any choice of coefficients
$\alpha_{l}(\bar{\bm{F}},\bm{\mu})$, all the $N\cdot L$ terms in Eq. 25 that
approximate the integral have to be approximated as accurately as possible.
Hence, the goal becomes to find a subset $Q(\ll\hat{Q})$ of integration points
with corresponding weights $\\{(\bm{X}_{q},w_{q})\\}_{q=1}^{Q}$ that
approximates Eq. 25 well, i.e.,

$\displaystyle\sum_{q=1}^{\hat{Q}}\hat{w}_{q}\left.\left[\frac{\partial\bm{\phi}_{i}}{\partial\bm{X}}:\bm{B}_{l}\right]\right|_{\hat{\bm{X}}_{q}}\approx\sum_{q=1}^{Q}w_{q}\left.\left[\frac{\partial\bm{\phi}_{i}}{\partial\bm{X}}:\bm{B}_{l}\right]\right|_{\bm{X}_{q}},\qquad
i=1,\dots,N,\ l=1,\dots,L.$ (26)

These $Q$ points and corresponding weights are found using a greedy algorithm,
the details of which can be found in [23] and are omitted here. The algorithm
is terminated when the mean squared error of all $N\cdot L$ terms is less than
a user-defined tolerance $\varepsilon$.

Compared to the original algorithm for a fixed geometry, as proposed in [23],
the only differences are that the weighted stress $\bm{W}$ is employed instead
of the stress $\bm{P}$ and that the parent domain $\Omega^{\rm{p}}$ is
considered instead of a fixed domain $\Omega$. With the ECM integration rule,
the hyper-reduced force vector and global stiffness matrix are computed as

$\displaystyle f^{\text{PODECM}}_{i}(\mathbf{a})$
$\displaystyle\coloneqq\sum_{q=1}^{Q}w_{q}\left.\left[\left(\frac{\partial\bm{\phi}_{i}}{\partial\bm{X}^{\rm{p}}}\bm{F}_{\bm{\mu}}^{-1}\right):\bm{P}\left(\bar{\bm{F}}+\left(\sum_{n=1}^{N}a_{n}\frac{\partial\bm{\phi}_{n}}{\partial\bm{X}^{\rm{p}}}\right)\bm{F}_{\bm{\mu}}^{-1}\right)\left|\det\bm{F}_{\bm{\mu}}\right|\right]\right|_{\bm{X}_{q}},$
(27) $\displaystyle K^{\text{PODECM}}_{ij}(\mathbf{a})$
$\displaystyle\coloneqq\sum_{q=1}^{Q}w_{q}\left.\left[\left(\frac{\partial\bm{\phi}_{i}}{\partial\bm{X}^{\rm{p}}}\bm{F}_{\bm{\mu}}^{-1}\right):\bm{A}\left(\bar{\bm{F}}+\left(\sum_{n=1}^{N}a_{n}\frac{\partial\bm{\phi}_{n}}{\partial\bm{X}^{\rm{p}}}\right)\bm{F}_{\bm{\mu}}^{-1}\right):\left(\frac{\partial\bm{\phi}_{j}}{\partial\bm{X}^{\rm{p}}}\bm{F}_{\bm{\mu}}^{-1}\right)\left|\det\bm{F}_{\bm{\mu}}\right|\right]\right|_{\bm{X}_{q}}.$
(28)

###  3.3 Effective Quantities

Once the new set of integration points and weights is found, the integrands of
Eqs. 27 and 28 only need to be evaluated at the points
$\\{\bm{X}_{q}\\}_{q=1}^{Q}$ during the solution of the reduced problem. This
also means that the stress and stiffness field are available at these points
only. To compute the effective quantities, the most straightforward method is
to use the integration rule obtained by ECM, i.e.,

$\displaystyle\begin{aligned}
\bar{\bm{P}}&=|\Omega^{\rm{p}}|^{-1}\int_{\Omega^{\rm{p}}}\bm{P}^{*\rm{p}}|\det\bm{F}_{\bm{\mu}}|d\bm{X}^{\rm{p}}\\\
&\approx|\Omega^{\rm{p}}|^{-1}\sum_{q=1}^{Q}w_{q}\left.\left(\bm{P}^{*\rm{p}}|\det\bm{F}_{\bm{\mu}}|\right)\right|_{\bm{X}_{q}}.\end{aligned}$
(29)

Since the stress field $\bm{P}^{*\rm{p}}$ is known at all integration points
$\\{\bm{X}_{q}\\}_{q=1}^{Q}$, the effective stress can be directly evaluated.
The method yields very accurate results in the examples considered below in
Section 4. However, it should be noted that there is currently no guarantee
that the integration rule found by ECM will generally be accurate for the
computation of the effective stress.

As discussed in Section 2.1.1, derivatives
$\dfrac{\partial\bm{w}^{*}}{\partial\bar{\bm{F}}}$ are needed to find the
effective stiffness $\bar{\bm{A}}$, see Eq. 12. For each component of
$\bar{\bm{F}}$, the linear tangent problem of Eq. 13 needs to be solved. By
employing the trial space of the fluctuation field for the auxiliary function
$\bm{q}_{kl}$, i.e.,

$\displaystyle\bm{q}_{kl}=\sum_{n=1}^{N}q_{n}\bm{\phi}_{n}(\bm{X}^{\rm{p}}),$
(30)

and the integration rule found by ECM, the following linear system results:

$\displaystyle\mathbf{K}^{*\rm{p}}\mathbf{q}=\mathbf{b},$ (31)

where $\mathbf{q}=[q_{1},\dots,q_{N}]^{T}$ is the column vector of unknowns to
be solved for and

$\displaystyle K^{*\rm{p}}_{ij}$
$\displaystyle=\sum_{q=1}^{Q}w_{q}\left.\left[\left(\frac{\partial\bm{\phi}_{i}}{\partial\bm{X}^{\rm{p}}}\bm{F}_{\bm{\mu}}^{-1}\right):\bm{A}^{*\rm{p}}:\left(\frac{\partial\bm{\phi}_{j}}{\partial\bm{X}^{\rm{p}}}\bm{F}_{\bm{\mu}}^{-1}\right)\left|\det\bm{F}_{\bm{\mu}}\right|\right]\right|_{\bm{X}_{q}},$
(32) $\displaystyle b_{i}$
$\displaystyle=-\left(\sum_{q=1}^{Q}w_{q}\left.\left[\left(\frac{\partial\bm{\phi}_{i}}{\partial\bm{X}^{\rm{p}}}\bm{F}_{\bm{\mu}}^{-1}\right):\bm{A}^{*\rm{p}}\left|\det\bm{F}_{\bm{\mu}}\right|\right]\right|_{\bm{X}_{q}}\right):\mathbb{E}_{kl}.$
(33)

Note that the matrix $\mathbf{K}^{*\rm{p}}\in\mathbb{R}^{N\times N}$ is
exactly the same as the global stiffness matrix $\mathbf{K}$ of Eq. 28
evaluated at the solution $\bm{w}^{*\rm{p}}$. After solving the tangent
problems, the effective stiffness $\bar{\bm{A}}$ can be computed (in index
notation) as

$\displaystyle\bar{A}_{ijkl}$
$\displaystyle=|\Omega^{\rm{p}}|^{-1}\sum_{q=1}^{Q}w_{q}\left.\left(\frac{\partial
P_{ij}^{*\rm{p}}}{\partial\bar{F}_{kl}}|\det\bm{F}_{\bm{\mu}}|\right)\right|_{\bm{X}_{q}},$
(34)

where

$\displaystyle\begin{aligned}
\frac{\partial\bm{P}^{*\rm{p}}}{\partial\bar{F}_{kl}}&=\bm{A}^{*\rm{p}}:\left(\mathbb{E}_{kl}+\left(\frac{\partial\bm{q}_{kl}}{\partial\bm{X}^{\rm{p}}}\bm{F}_{\bm{\mu}}^{-1}\right)\right).\end{aligned}$
(35)

The derivation of the effective sensitivities (recall Eq. 14) can be performed
analogously to the effective stiffness. However, due to the derivatives of the
inverse and determinant of the transformation gradient $\bm{F}_{\bm{\mu}}$,
such derivations become rather technical and hence are omitted.

###  3.4 Auxiliary Problem For Geometrical Transformation

Thus far, the geometrical transformation
$\bm{\Phi}_{\bm{\mu}}:\Omega^{\rm{p}}\rightarrow\Omega^{\bm{\mu}}$ has been
assumed to be known and has not been discussed in more detail. However, such
transformations are in general not known analytically and have to be found
numerically by using, for example, radial basis functions, see, e.g., [30,
31], or mesh-based methods, see, e.g., [32, 29]. For each of those methods, an
auxiliary problem arises which needs to be solved. In order to rapidly solve
the surrogate model for a wide range of different geometries, it must
therefore also be ensured that the auxiliary problem can be solved rapidly. In
this work, the method in [29] is employed, in which the auxiliary problem is
formulated as a linear elasticity problem by defining
$\bm{\Phi}_{\bm{\mu}}(\bm{X}^{\rm{p}})=\bm{X}^{\rm{p}}+\bm{d}(\bm{X}^{\rm{p}})$,
with $\bm{d}$ the transformation displacement obtained from

$\displaystyle\text{Div}\left(\mathbb{C}^{\text{aux}}:\frac{1}{2}\left(\frac{\partial\bm{d}}{\partial\bm{X}^{\rm{p}}}+\left(\frac{\partial\bm{d}}{\partial\bm{X}^{\rm{p}}}\right)^{T}\right)\right)$
$\displaystyle=\bm{0}$ $\displaystyle\text{ in }\Omega^{\rm{p}}.$ (36)

In the above equation, $\mathbb{C}^{\text{aux}}$ is the fourth-order
elasticity tensor, fully specified by the Young’s modulus $E^{\text{aux}}$ and
Poisson’s ratio $\nu^{\text{aux}}$. The boundary conditions for this PDE are
problem-dependent and are specified by the geometrical parameters $\bm{\mu}$.
For the RVE problem, the outer boundaries are fixed ($\bm{d}=\bm{0}$), while
$\bm{d}$ is prescribed on parts of the interior that are parameterized by
$\bm{\mu}$. In [29] the effect of the choice of $E^{\text{aux}}$ and
$\nu^{\text{aux}}$ was studied and it was demonstrated empirically that the
choice only has a minor effect on the final approximation quality. Hence, in
all numerical examples considered in this work, a Young’s modulus of
$E^{\text{aux}}=1$ and Poisson’s ratio $\nu^{\text{aux}}=0.25$ is assumed. The
auxiliary problem can then be significantly accelerated with the RBM in
combination with a DEIM (e.g., [29, 31]), resulting in

$\displaystyle\hat{\mathbf{A}}\hat{\mathbf{d}}=\hat{\mathbf{b}}(\bm{\mu}),$
(37)

where $\hat{\mathbf{A}}\in\mathbb{R}^{N_{p}\times N_{p}}$ is the reduced
system matrix, $\hat{\mathbf{d}}\in\mathbb{R}^{N_{p}}$ is the reduced
transformation displacement, $\hat{\mathbf{b}}(\bm{\mu})\in\mathbb{R}^{N_{p}}$
is the reduced forcing vector and $N_{p}$ is the number of geometrical
parameters. Since $N_{p}$ is usually small, Eq. 37 can be rapidly solved. From
$\hat{\mathbf{d}}$, the transformation gradient $\bm{F}_{\bm{\mu}}$, its
inverse $\bm{F}_{\bm{\mu}}^{-1}$, and its determinant $\det\bm{F}_{\bm{\mu}}$
can be computed. Moreover, expressions for the derivative of the inverse and
determinant of the transformation gradient $\bm{F}_{\bm{\mu}}$ can be derived,
which are needed for computing the sensitivities with respect to $\bm{\mu}$,
see Eq. 14. For more information on the auxiliary problem and its reduction,
the reader is referred to [29].

###  3.5 Summary

For convenience, the offline-online decomposition for constructing and solving
the surrogate model is summarized in Algorithm 1.

1:

2:Define a parent domain $\Omega^{\rm{p}}$ and its finite element
discretization.

3:Generate parameter samples
$\\{\bar{\bm{F}}^{i},\bm{\mu}^{i}\\}_{i=1}^{N_{s}}$ from a random
distribution.

4:For each different set of geometrical parameters $\bm{\mu}^{i}$, solve the
auxiliary problem in Eq. 36 to obtain the transformation map
$\bm{\Phi}_{\bm{\mu}^{i}}$.

5:Compute $\bm{F}_{\bm{\mu}^{i}}^{-1}$ and $\det\bm{F}_{\bm{\mu}^{i}}$ for
each parameter sample $\bm{\mu}^{i}$, then run full simulations (Eqs. 5 and 6)
for $\bar{\bm{F}}^{i}$ and collect fluctuation displacement and weighted
stress snapshots.

6:Compute POD for the fluctuation displacement and weighted stress, cf. Eqs.
15 and 23.

7:Run ECM algorithm and find integration points and weights, cf. Eq. 26.

8:Assemble the reduced system matrix and forcing vector for the auxiliary
problem in Eq. 37 by applying POD and DEIM. Details are provided in [29].

9:

1:Given a new parameter set $(\bar{\bm{F}}^{*},\bm{\mu}^{*})$, solve reduced
auxiliary problem Eq. 37 and compute $\bm{F}_{\bm{\mu}^{*}}^{-1}$ and
$\det\bm{F}_{\bm{\mu}^{*}}$.

2:Solve reduced problem for $\bar{\bm{F}}^{*}$ with Eqs. 27 and 28.

3:Compute effective stress using Eq. 29.

4:Solve the linear problem Eq. 31 for each component of $\bar{\bm{F}}^{*}$.

5:Compute components of the effective stiffness with Eq. 34.

Algorithm 1 Offline-online decomposition of the proposed PODECM framework with
microstructures parameterized with external loading $\bar{\bm{F}}$ and
geometrical features $\bm{\mu}$.

##  4 Example Problems

The proposed framework, referred to as PODECM, is first tested on a non-linear
composite microstructure under various loading conditions and analyzed in
depth regarding its capabilities and accuracy. The RVE consists of an elasto-
plastic matrix with stiff inclusions of variable size and is considered under
non-monotonic loading. The surrogate model is analyzed in terms of the number
of basis functions of the fluctuation displacement field $N$, number of basis
functions of the weighted stress $L$ and the ECM integration error tolerance
$\varepsilon$. Subsequently, a two-scale problem involving a porous
microstructure under non-monotonic loading conditions and varying porosities
is studied to illustrate the accuracy and speed-up of PODECM in a two-scale
setting.

All experiments are defined in two dimensions under plane strain conditions.
The RVEs are assumed to be of size $[0,1]^{2}$ and all quantities are assumed
to be normalized and hence dimensionless. Since the macroscopic deformation
gradient $\bar{\bm{F}}$ can always be decomposed into a rotation
$\bar{\bm{R}}$ and a symmetric stretch tensor $\bar{\bm{U}}$ with a polar
decomposition, i.e., $\bar{\bm{F}}=\bar{\bm{R}}\bar{\bm{U}}$, it is sufficient
to generate training data for the stretch tensor $\bar{\bm{U}}$, having only 3
independent components (6 in 3D).

To measure the quality of the approximation, the following error measures to
compare the full FE simulations against PODECM solutions are defined:

  1. 1.

Error of effective stress

$\displaystyle\epsilon_{\bar{\bm{P}}}=\frac{\sum_{k=1}^{K}||\bar{\bm{P}}^{\text{PODECM}}(\bar{\bm{U}}^{k})-\bar{\bm{P}}^{\text{FE}}(\bar{\bm{U}}^{k})||_{F}}{\sum_{k=1}^{K}||\bar{\bm{P}}^{\text{FE}}(\bar{\bm{U}}^{k})||_{F}},$
(38)

where $\bar{\bm{P}}^{\text{PODECM}}(\bar{\bm{U}}^{k})$ and
$\bar{\bm{P}}^{\text{FE}}(\bar{\bm{U}}^{k})$ denote the effective stress
obtained with PODECM and FE for $\bar{\bm{U}}^{k}$, $||\bullet||_{F}$ denotes
the Frobenius norm, $K$ is the total number of loading steps and
$\bar{\bm{U}}^{k}$ is the applied external load at load step $k$.

  2. 2.

Error of fluctuation field

$\displaystyle\epsilon_{\bm{w}}=\frac{\sum_{k=1}^{K}||\bm{w}^{\text{PODECM}}(\bar{\bm{U}}^{k})-\bm{w}^{\text{FE}}(\bar{\bm{U}}^{k})||_{\mathcal{V}}}{\sum_{k=1}^{K}||\bm{w}^{\text{FE}}(\bar{\bm{U}}^{k})||_{\mathcal{V}}},$
(39)

where $\bm{w}^{\text{PODECM}}$ and $\bm{w}^{\text{FE}}$ denote the fluctuation
displacement field obtained with PODECM and FE, and
$||\bullet||_{\mathcal{V}}^{2}=(\bullet,\bullet)_{\mathcal{V}}$, c.f., Eq. 4.
Recall that the integral in Eq. 4 is defined over the parameterized domain
$\Omega^{\bm{\mu}}$.

###  4.1 Elasto-Plastic Composite RVE With Random Stiff Inclusions

####  4.1.1 Problem Description

The considered RVE in this example consists of two phases, an elasto-plastic
matrix and stiff elastic inclusions. The geometry of the parent domain is
shown in Fig. 2(a), where the volume fraction of the inclusions is $23.4\%$.
For the geometrical parameterization, one geometrical parameter
$\bm{\mu}=\\{\zeta\\}$ that scales the size of the inclusions uniformly (and
is proportional to the volume fraction of the inclusions) is introduced, see
Figs. 2(b) and 2(c) showing two example domains for distinct values of
$\zeta$. The simulation mesh is depicted in Fig. 2(d), where six-noded
quadratic triangular elements are used in conjunction with three quadrature
points per element. In total, the mesh has 62194 degrees of freedom, 15450
triangular elements and 46350 quadrature points.

(a) Parent domain $\Omega^{\rm{p}}$ with $\zeta=1$

(b) Domain $\Omega^{\bm{\mu}_{1}}$ with $\zeta=1.2$

(c) Domain $\Omega^{\bm{\mu}_{2}}$ with $\zeta=0.5$

(d) Simulation mesh

Figure 2: Parent with two parameterized domains and simulation mesh. (a) The
parent domain consists of a matrix material (blue) with 23 random elliptical
inclusions (orange). The problem has one geometrical parameter $\zeta$ that
scales all ellipses uniformly ($\zeta=1$ for parent domain). (b) A
parameterized domain for $\zeta=1.2$ and (c) for $\zeta=0.5$. (d) The
considered mesh consists of six-noded triangular elements and contains in
total 62194 degrees of freedom, 15450 triangular elements and 46350 quadrature
points.

For the constitutive model of both matrix and inclusion the small-strain
$J_{2}$-plasticity model with linear isotropic hardening is chosen and
extended to large strains with the method presented in Cuitino and Ortiz [33].
The small-strain model obeys:

$\displaystyle\bm{\sigma}$
$\displaystyle=\mathbb{D}:(\bm{\epsilon}-\bm{\epsilon}_{\text{pl}}),$ (40)
$\displaystyle f^{\text{yield}}$
$\displaystyle=||\bm{\sigma}||_{\text{mises}}-(\sigma_{y0}+H\xi),$ (41)
$\displaystyle\bm{r}$ $\displaystyle=\frac{\partial
f}{\partial\bm{\sigma}}=\sqrt{\frac{3}{2}}\frac{\operatorname{Dev}(\bm{\sigma})}{\sqrt{\operatorname{Dev}(\bm{\sigma}):\operatorname{Dev}(\bm{\sigma})}},$
(42) $\displaystyle\dot{\xi}$ $\displaystyle=\gamma,$ (43)
$\displaystyle\dot{\bm{\epsilon}}_{\text{pl}}$ $\displaystyle=\gamma\bm{r},$
(44) $\displaystyle\gamma$ $\displaystyle\geq 0,\ f^{\text{yield}}\geq 0,\
\gamma f^{\text{yield}}=0,$ (45)

where $\bm{\epsilon}$ is the small-strain tensor, $\bm{\epsilon}_{\text{pl}}$
the plastic strain tensor, $\mathbb{D}$ is the fourth-order elasticity tensor
that can be fully specified by Young’s modulus $E$ and Poisson’s ratio $\nu$,
and $\bm{\sigma}$ is the corresponding stress tensor; $f^{\text{yield}}$
defines the yield surface, $\operatorname{Dev}(\bullet)$ takes the deviatoric
part of a tensor $(\bullet)$,
$||\bullet||_{\text{mises}}=\sqrt{\dfrac{3}{2}\operatorname{Dev}(\bullet):\operatorname{Dev}(\bullet)}$
computes the von Mises stress, $H$ is the hardening constant, $\sigma_{y0}$
yield stress, $\bm{r}$ the plastic flow direction, $\xi$ the equivalent
plastic strain that defines the isotropic hardening of the yield surface, and
$\gamma$ is the consistency parameter. For more information on the material
model, see Simo and Hughes [34].

Following the procedure of [33], by employing a multiplicative split of the
deformation gradient $\bm{F}=\bm{F}^{\text{el}}\bm{F}^{\text{pl}}$, the
elastic logarithmic strain can be defined as

$\displaystyle\bm{C}^{\text{el}}_{\text{log}}\coloneqq\ln{\bm{C}^{\text{el}}}=\ln((\bm{F}^{\text{el}})^{T}\bm{F}^{\text{el}}).$
(46)

By interpreting the elastic logarithmic strain as the small-strain tensor, the
small-strain constitutive model defined in Eqs. 40, 41, 42, 43, 44 and 45 is
used to compute the stress tensor $\hat{\bm{S}}$ on the intermediate
configuration,

$\displaystyle\hat{\bm{S}}\coloneqq
2\bm{\sigma}(\bm{C}^{\text{el}}_{\text{log}}):\dfrac{\partial{\bm{C}^{\text{el}}_{\text{log}}}}{\partial{\bm{C}^{\text{el}}}},$
(47)

while the 1PK stress is recovered from $\hat{\bm{S}}$ as

$\displaystyle\bm{P}=(\bm{F}^{\text{el}})^{-1}\hat{\bm{S}}\
(\bm{F}^{\rm{pl}})^{-T}.$ (48)

Instead of evolving the plastic strain with Eq. 44, the plastic deformation
gradient $\bm{F}^{\rm{pl}}$ is evolved according to

$\displaystyle\dot{\bm{F}}^{\rm{pl}}=\exp(\gamma\bm{r})\bm{F}^{\rm{pl}}.$ (49)

For the matrix, the following material parameters are selected: $E=10$,
$\nu=0.3$, $\sigma_{y0}=0.2$ and $H=5$. For the inclusions, $E=100$ and
$\nu=0.3$ is selected, corresponding to a stiffness contrast ratio of 10.
Since no plastic deformation is assumed for the inclusions, its yield stress
is set to a large value such that yielding never occurs.

Three loading $\bar{U}_{xx}$, $\bar{U}_{xy}$, $\bar{U}_{yy}$ and one
geometrical $\zeta$ parameters are considered with bounds $\zeta\in[0.5,1.2]$,
$\bar{U}_{xx}\in[0.9,1.1]$, $\bar{U}_{yy}\in[0.9,1.1]$ and
$\bar{U}_{xy}\in[-0.1,0.1]$. Through $\zeta$, the volume fraction of the
inclusions is varied from 5.85% to 33.7%. The macroscopic stretch tensor is
applied in $K=40$ load steps. For the first 10 load steps the applied stretch
increases linearly from $\bm{0}$ to $\bar{\bm{U}}$. In the next 20 steps, the
applied loading decreases linearly from $\bar{\bm{U}}$ to $-\bar{\bm{U}}$,
such that at $k=20$ the applied stretch is $\bm{0}$. Finally, in the last 10
steps, the stretch increases linearly from $-\bar{\bm{U}}$ to $\bm{0}$. Even
though only macroscopic strains of up to 10% are applied, local strains
reached values up to 83%. Evolution of the effective von Mises stress for an
example with $\zeta=1.010$, $\bar{U}_{xx}=1.1$, $\bar{U}_{yy}=1.0$,
$\bar{U}_{xy}=0.0$ is shown in Fig. 3, including the local von Mises stress
fields at steps $k=\\{10,20,30,40\\}$.

Figure 3: Macroscopic von Mises stress $\bar{\bm{\sigma}}_{\text{mises}}$
plotted over $k$ for $\zeta=1.010$, $\bar{U}_{xx}=1.1$, $\bar{U}_{yy}=1.0$,
$\bar{U}_{xy}=0.0$. At $k=\\{10,20,30,40\\}$, microstrucural von Mises stress
fields $\bm{\sigma}_{\text{mises}}$ are shown. The von Mises stress is non-
zero at $k=\\{20,40\\}$ due to residual plastic deformation.

####  4.1.2 Results

Figure 4: The left column shows the number of quadrature points $Q$ obtained
from ECM for different choices of number of basis functions of fluctuation
field $N$, number of basis functions of weighted stress $L$, and ECM
integration error $\varepsilon$. The middle and right columns show the average
errors of the effective stress and the fluctuation field when tested on the
testing data for different choices of $N$, $L$ and $\varepsilon$. The top row
assumes a fixed $\varepsilon=0.01$, while the bottom row assumes a fixed
$L=15$.

In total 20 samples are generated from a Sobol sequence to train PODECM
whereas 100 testing samples are generated from a uniform distribution to test
PODECM. Each sample consists of 40 snapshots for each load step.

The accuracy and speed-up of PODECM depends on the number of basis functions
used for the fluctuation displacement field $N$ and the number of quadrature
points $Q$. While $N$ is typically chosen directly, $Q$ depends on the choice
of the number of basis functions used for the weighted stress $L$ and the ECM
integration error $\varepsilon$.

To study the influence of $L$ on the resulting number of quadrature points $Q$
and mean errors in effective stress and fluctuation field on the testing
dataset, several combinations of $N$ and $L$ for a fixed $\varepsilon=10^{-2}$
are tested, with resulting errors shown in the top row of Fig. 4. The
projection error (for $N$ basis functions and using full integration) is shown
as well. It can be clearly seen that the number of quadrature points $Q$
increases drastically with increasing $N$ and $L$, as more information needs
to be integrated accurately. For the mean errors, a higher $L$ leads to better
results on average, although we observe that errors fluctuate significantly,
and for some values of $N$ a worse approximation is obtained with a higher
$L$. This occurs since the ECM algorithm is a greedy algorithm, meaning that
it does not necessarily find an optimal set of integration points. When more
basis functions are included into the algorithm, a completely different set of
points may be found that finally leads to a worse approximation. It can
furthermore be observed that the gap between the projection error and the
PODECM solution grows larger for increasing $N$. This is because the basis
functions typically become more oscillatory and difficult to approximate with
higher $N$, see, e.g., [35, 36], and thus require significantly more
quadrature points for a good approximation. It is interesting that the gap for
the errors in the fluctuation field are smaller than the ones in the effective
stress. This happens because the ECM integration rule used to compute the
effective stress introduces an additional approximation and hence also a
source of error.

Several combinations of $N$ and $\varepsilon$ for a fixed $L=15$ are next
tested to study the influence of $\varepsilon$ on the number of quadrature
points and approximation errors. Obtained results are shown in the bottom row
of Fig. 4. Similarly to the previous analysis, a lower $\varepsilon$ leads to
more quadrature points $Q$ and a lower mean error in the effective stress and
fluctuation field on average, as the integrals are approximated more
accurately. Interestingly, lowering the tolerance from 0.01 to 0.001 does not
significantly improve the approximation quality, even though substantially
more quadrature points are included, meaning that the errors can be attributed
to the higher modes of the weighted stress (the additional quadrature points
barely contain any information). Therefore, choosing a tolerance smaller than
$\varepsilon=0.01$ leads to no improvement.

From Fig. 4 we further observe that the errors of the fluctuation field are
considerably higher (order of magnitude) than the errors of the effective
stresses. This results from the fact that the POD basis functions aim to
minimize the $H_{1}(\Omega^{\rm{p}})$ error, and thus approximate the field
accurately on average rather than locally, suggesting a favorable
approximation for averaged quantities such as the effective stress. In [37],
the authors showed a similar result for FFT: the effective stresses converge
with an order of $h$, while the local fields converge with $h^{1/2}$, where
$h$ is the used voxel edge-length.

To conclude, the more basis functions $N$ and $L$ are used and the lower the
integration error $\varepsilon$ is chosen, the more accurate the final result
is. However, at the same time the surrogate model grows in size and the speed-
up decreases. A user must thus make a compromise between accuracy and cost. In
our experience, the speed-up correlates nearly linearly with $Q$, i.e., if the
number of quadrature points is reduced by a factor of 100, this results in a
speed-up of roughly 100. In contrast, the number of basis functions $N$ only
plays a minor role for the speed-up. For this example, use of $N=20$, $L=20$
and $\varepsilon=0.01$ lead to a reduction in the number of degrees of freedom
from 62194 to 20 and quadrature points from 46350 to 212, suggesting a speed-
up on the order of roughly 200.

###  4.2 Two-Scale Compression With Porous Microstructure

####  4.2.1 Problem Description

In the second example, the macroscopic structure, depicted in Fig. 5 together
with the employed simulation mesh, is compressed under an external loading
$T(x)$. Here, we assume $H=1$, $W=2$,
$T(x)=\bar{T}\left(1-\left(\dfrac{2x}{W}-1\right)^{2}\right)$, $x\in[0,W]$,
with $\bar{T}$ the magnitude of the applied load. The simulation mesh has 1322
degrees of freedom, 200 8-noded quadrilateral elements, and in total 800
quadrature points. The structure is assumed to have a porous microstructure,
modelled by the parameterized RVE, shown in Fig. 6, and the same material
model as the matrix material in the previous example Eqs. 40, 41, 42, 43, 44,
45, 46, 47, 48 and 49. Such microstructures (with circular holes, i.e., $a=b$)
have been considered in several works, see, e.g., [38, 39], due to their
auxetic behavior under compression, i.e., negative Poisson’s ratio. During
compression, the center part of the material starts to rotate, thus pulling
the material from the sides inwards. In our work, we define two independent
parameters, namely the volume fraction of the voids, $v_{\text{void}}\coloneqq
4\pi ab$, and the ratio of the semi-major axis $b$ to semi-minor axis $a$ of
each hole, $\kappa\coloneqq b/a$. The semi-minor axis $a$ and semi-major axis
$b$ depend on $v_{\text{void}}$ and $\kappa$ as

$\displaystyle a(v_{\text{void}},\kappa)$
$\displaystyle=\sqrt{v_{\text{void}}/(4\pi\kappa)},$ (50) $\displaystyle
b(v_{\text{void}},\kappa)$ $\displaystyle=\kappa\cdot
a(v_{\text{void}},\kappa).$ (51)

(a) Macroscopic structure

(b) Macroscopic simulation mesh

Figure 5: Geometry (a) and mesh (b) of the considered macroscopic structure.
The body is fixed on the bottom and an external compression force $T$ is
applied on the top. The mesh consists of 1322 degrees of freedom, 200 8-noded
quadrilateral elements and 4 quadrature points per element.

(a) RVE

(b) Simulation mesh

Figure 6: Geometry (a) and mesh (b) of the porous RVE. The elliptical holes
have semi-minor axis $a$ and semi-major axis $b$, and are parameterized by the
volume fraction of the pores $v_{\text{void}}$ and the $\kappa=b/a$. The
employed simulation mesh has 21042 degrees of freedom, 4964 6-noded triangular
elements and in total 14892 quadrature points. The parent domain corresponds
to $\kappa=1.25$ and $v_{\text{void}}=0.45$.

(a) Effective Poisson’s ratio $\nu^{\text{eff}}$

(b) Effective Young’s modulus $E^{\text{eff}}$

Figure 7: Initial effective Poisson’s ratio (a) and Young’s modulus (b) of RVE
for different values of $v_{\text{void}}$ and $\kappa$.

Depending on the values of the parameters, the resulting effective properties
change significantly. To illustrate this, linear analyses of this RVE for
different parameters have been carried out, similarly to [40], where a small
compression in the $y$-direction with $\Delta u_{y}=0.001$ has been applied,
while allowing the RVE to contract freely in the $x$-direction. With the
resulting displacement in the $x$-direction, $\Delta u_{x}$, the Poisson’s
ratio in the initial state can be estimated as

$\displaystyle\nu^{\text{eff}}=-\frac{\Delta u_{x}}{\Delta u_{y}}.$ (52)

Similarly, the initial Young’s modulus is estimated as

$\displaystyle E^{\text{eff}}=\frac{\bar{P}_{yy}}{\Delta u_{y}},$ (53)

where $\bar{P}_{yy}$ is the $yy$-component of the effective stress. For
parameter ranges $v_{\text{void}}\in[0.4,0.5]$ and $\kappa\in[1.01,1.5]$, the
estimated Poisson’s ratio and Young’s modulus are plotted in Fig. 7. It can be
observed that removing material (by increasing $v_{\text{void}}$ while keeping
$\kappa$ fixed) or increasing $\kappa$ while keeping $v_{\text{void}}$ fixed
both lead to a softer response with lower Young’s modulus. While the Poisson’s
ratio is barely affected by $v_{\text{void}}$ for values of $\kappa$ close to
1, the effect becomes apparent for larger values of $\kappa$. In particular,
for $\kappa\geq 1.4$, the Poisson’s ratio changes from a positive value to a
negative one. Therefore, by tuning $v_{\text{void}}$ and $\kappa$, the RVE
behavior can be significantly modified.

The parameters $v_{\text{void}}$ and $\kappa$ are chosen to vary smoothly
through the macrostructural domain as

$\displaystyle\kappa(x,y)$ $\displaystyle=1.5-(1.5-1.01)y,$ (54)
$\displaystyle v_{\text{void}}(x,y)$ $\displaystyle=0.4+(0.5-0.4)(1-x)^{2}.$
(55)

The parent mesh is selected with $\kappa=1.25$ and $v_{\text{void}}=0.45$. The
external loading is applied in $K=50$ load steps. In the first 25 load steps,
the applied load $\bar{T}$ increases linearly from 0 to 0.2. In the next 25
steps $\bar{T}$ is decreased linearly from 0.2 to 0.

####  4.2.2 Results

Several two-scale simulations with different PODECM surrogate models are run
and compared to the full reference FE2 solution. To compare the accuracy of
the surrogate solutions, the compliance
$C\coloneqq\int_{\Gamma}T(x)u_{y}(x)dx$, where $\Gamma$ denotes the top
horizontal edge of the macrostructure and $u_{y}$ its vertical displacement,
is computed at every load step. The compliance is an important quantity, often
employed in optimization problems. Subsequently, the relative error in
compliance $\epsilon_{C}$ and the relative error averaged over all load steps
$\bar{\epsilon}_{C}$ are defined as

$\displaystyle\epsilon_{C,k}\coloneqq\frac{|C_{k}-C_{k}^{\text{FE}^{2}}|}{|C_{k}^{\text{FE}^{2}}|},\qquad\bar{\epsilon}_{C}\coloneqq\frac{1}{K}\sum_{n=1}^{K}\epsilon_{C,k},$
(56)

where the subscript $k$ denotes the $k$-th load step and $C^{\text{FE}^{2}}$
is the compliance computed with the full solution.

The tested PODECM RVE models are generated for different numbers of training
samples $N_{\text{train}}$ and number of basis functions $N$. The training
data is sampled from $\bar{U}_{xx}\in[0.85,1]$, $\bar{U}_{yy}\in[0.85,1]$,
$\bar{U}_{xy}\in[-0.15,0.15]$, $v_{\text{void}}\in[0.4,0.5]$ and
$\kappa\in[1.01,1.5]$ with a Sobol sequence. Each sample consists of 50 load
steps, where the macroscopic stretch tensor is linearly increased for 25 load
steps and then linearly decreased to 0 in 25 steps. For the ECM algorithm, the
number of basis functions for the weighted stress $L$ and the integration
error $\varepsilon$ are assumed as fixed with $L=20$ and $\varepsilon=0.01$
for all models. The exact settings for each surrogate model are summarized in
Table 1, alongside the averaged relative error in compliance
$\bar{\epsilon}_{C}$, and the run time and speed-up in comparison to the full
FE2 solution. For comparison, the total number of degrees of freedom and
quadrature points of the full FE model of the RVE are also provided. As can be
seen, for all surrogate models the number of quadrature points are reduced by
a factor of up to 100 which results in high speed-ups up to 100 times, while
errors are below 5% for all models. By increasing $N$ from 10 to 50 for rom_1
to rom_5, the error decreases from 4.74% to 1.54%, whereas the speed-up
reduces from 92 to 26. Including more training samples with the same number of
basis functions $N=50$ (from $N_{\text{train}}=20$ for rom_5 to
$N_{\text{train}}=100$ for rom_6) improves the error from 1.54% to 0.65% while
the speed-up remains roughly the same. This means that by increasing the
sample size, the first 50 basis functions contain more general information
that result in a better approximation.

Table 1: Summary of results for full FE2 and different PODECM surrogate models. Reduced order models are generated for different numbers of training samples $N_{\text{train}}$ and basis functions of the displacement $N$. The number of quadrature points $Q$ follow from the ECM algorithm with a fixed number of basis functions of the weighted stress field $L=20$ and an integration error of $\varepsilon=0.01$. All reduced order models achieve errors less than 5% with speed-ups up to 100 times. By generating more training data and maintaining the same $N$ (rom_5 and rom_6), better results are achieved. | $N_{\text{train}}$ | $N$ | $Q$ | $\bar{\epsilon}_{C}$ | run time111All computations are executed using 20 cores of an Intel Platinum 8260. | speed-up
---|---|---|---|---|---|---
full | - | 21042 | 14892 | - | 12573s | -
rom_1 | 20 | 10 | 132 | 4.74% | 136s | 92.45
rom_2 | 20 | 20 | 259 | 4.66% | 205s | 61.33
rom_3 | 20 | 30 | 372 | 2.28% | 284s | 44.27
rom_4 | 20 | 40 | 490 | 2.19% | 385s | 32.66
rom_5 | 20 | 50 | 595 | 1.54% | 488s | 25.76
rom_6 | 100 | 50 | 577 | 0.65% | 496s | 25.34

In Fig. 8 the force displacement curve is shown for the FE2 and a few selected
surrogate solutions. The displacement is defined as the vertical displacement
at the mid point of the top edge (which is also the maximal displacement). It
can be observed that all surrogate models underpredict the displacement,
indicating that the surrogate models overpredict the stiffness of the
macrostructure.

Figure 8: Force-displacement curve of two-scale simulations. The displacement
$\tilde{u}$ is defined as the vertical displacement of the mid point on the
top edge, cf. Fig. 5(a). The structure starts deforming plastically for
$\bar{T}>0.07$ and a residual displacement of roughly $\tilde{u}=0.04$ remains
after unloading. All surrogate models achieve accurate results, although they
generally predict a slightly stiffer response than the full FE2 solution.

The relative error in compliance is plotted over the load steps $k$ in Fig. 9.
For all models, the error slowly increases over $k$. The reason for this
behavior is that all training samples are generated for simple loading cases,
where the macroscopic stretch tensor is linearly varied in only one direction
throughout the entire simulation. In the macroscopic simulation, the
macroscopic stretch tensor for one integration point generally does not evolve
along one direction, but changes its direction continuously, leading to highly
complicated deformation paths and histories that are not included in the
training data. To tackle this problem, random loading paths during training
could be used, as performed in, e.g., [41, 42], to generate a more general
surrogate model. Solely increasing the number of samples from 20 to 100 (rom_5
to rom_6) also decreases the observed errors to less than 1% for all load
steps. This shows that PODECM generalizes well to loading paths that are not
part of the training data.

Figure 9: Relative error $\epsilon_{C}$ in compliance over load step $k$. All
surrogate models (specified in Table 1) begin with a lower error which slowly
grows with $k$. By increasing the sample size of the training data (comparing
rom_5 and rom_6), the prediction becomes more accurate for the same number of
basis functions; the errors for all load steps for rom_6 are below 1%.

##  5 Conclusions

In this work, we developed a reduced order model, termed PODECM, by combining
the proper orthogonal decomposition (POD), the empirical cubature method (ECM)
and a geometrical transformation method. This method is designed to accelerate
microscopic computations within two-scale simulations, parameterized by
material and geometrical parameters, hence enabling PODECM to be used for two-
scale optimization problems. We showed how to compute the effective stress and
the corresponding consistent effective stiffness for this model, as required
by the macroscopic solver, and how to obtain the effective sensitivities with
respect to geometrical parameters. The framework was first tested on a single-
scale problem involving an RVE of a composite microstructure that consisted of
a soft elasto-plastic matrix with stiff inclusions of variable size,
controlled by a single geometrical parameter. With PODECM, the number of
degrees of freedom and integration points was reduced to a fraction of the
full FE model while maintaining a high accuracy in effective stress. The
performance of PODECM was further evaluated for a two-scale simulation, in
which a porous microstructure, characterized by two geometrical parameters,
was considered. Both geometrical parameters were varied throughout the
macrostructure, and depending on their values, the effective Poisson’s ratio
changed from positive to negative. For this example, different PODECM models
were constructed with good accuracies (of errors less than 1%), while
achieving speed-ups up to 100. Even though highly accurate and fast solutions
were obtained with the proposed method, several open questions prevail. For
instance, optimality and accuracy of the ECM integration rule cannot be
ensured and thus needs to be further analyzed and understood.

An attractive feature of the proposed framework is its versatility and
generality. Moreover, since the underlying microscopic PDE still needs to be
solved, only very little training data is required to construct a good
approximation of the full model with no special treatment of history-dependent
material behavior, which makes this framework viable for two-scale shape
optimization problems.

## Data availability

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

## Acknowledgements

This result is part of a project that has received funding from the European
Research Council (ERC) under the European Union’s Horizon 2020 Research and
Innovation Programme (Grant Agreement No. 818473). The authors would in
addition like to thank Martin Horák from Czech Technical University in Prague
for his help with the implementation of the large-strain $J_{2}$-plasticity
model.

## References

  * [1] M. G. D. Geers, V. G. Kouznetsova, W. A. M. Brekelmans, Multi-scale computational homogenization: Trends and challenges, Journal of Computational and Applied Mathematics 234 (7) (2010) 2175–2182.  doi:10.1016/j.cam.2009.08.077.
URL https://linkinghub.elsevier.com/retrieve/pii/S0377042709005536

  * [2] V. G. Kouznetsova, W. A. M. Brekelmans, F. P. T. Baaijens, An approach to micro-macro modeling of heterogeneous materials, Computational Mechanics 27 (1) (2001) 37–48.  doi:10.1007/s004660000212.
URL http://link.springer.com/10.1007/s004660000212

  * [3] C. Farhat, F.-X. Roux, A method of finite element tearing and interconnecting and its parallel solution algorithm, International Journal for Numerical Methods in Engineering 32 (6) (1991) 1205–1227.  doi:10.1002/nme.1620320604.
URL https://onlinelibrary.wiley.com/doi/10.1002/nme.1620320604

  * [4] P. Ladevèze, J.-C. Passieux, D. Néron, The LATIN multiscale computational method and the Proper Generalized Decomposition, Computer Methods in Applied Mechanics and Engineering 199 (21-22) (2010) 1287–1296.  doi:10.1016/j.cma.2009.06.023.
URL https://linkinghub.elsevier.com/retrieve/pii/S0045782509002643

  * [5] H. Moulinec, P. Suquet, A numerical method for computing the overall response of nonlinear composites with complex microstructure, Computer Methods in Applied Mechanics and Engineering 157 (1-2) (1998) 69–94.  doi:10.1016/S0045-7825(97)00218-1.
URL https://linkinghub.elsevier.com/retrieve/pii/S0045782597002181

  * [6] N. Mishra, J. Vondřejc, J. Zeman, A comparative study on low-memory iterative solvers for FFT-based homogenization of periodic media, Journal of Computational Physics 321 (2016) 151–168.  doi:10.1016/j.jcp.2016.05.041.
URL https://linkinghub.elsevier.com/retrieve/pii/S0021999116301863

  * [7] G. J. Dvorak, Transformation field analysis of inelastic composite materials, Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences 437 (1900) (1992) 311–327.  doi:10.1098/rspa.1992.0063.
URL https://royalsocietypublishing.org/doi/10.1098/rspa.1992.0063

  * [8] J. C. Michel, P. Suquet, Nonuniform transformation field analysis, International Journal of Solids and Structures 40 (25) (2003) 6937–6955.  doi:10.1016/S0020-7683(03)00346-9.
URL https://linkinghub.elsevier.com/retrieve/pii/S0020768303003469

  * [9] Z. Liu, M. A. Bessa, W. K. Liu, Self-consistent clustering analysis: An efficient multi-scale scheme for inelastic heterogeneous materials, Computer Methods in Applied Mechanics and Engineering 306 (2016) 319–341.  doi:10.1016/j.cma.2016.04.004.
URL http://dx.doi.org/10.1016/j.cma.2016.04.004

  * [10] C. Yu, O. L. Kafka, W. K. Liu, Self-consistent clustering analysis for multiscale modeling at finite strains, Computer Methods in Applied Mechanics and Engineering 349 (2019) 339–359.  doi:10.1016/j.cma.2019.02.027.
URL https://linkinghub.elsevier.com/retrieve/pii/S0045782519300982

  * [11] F. Feyel, Multiscale FE2 elastoviscoplastic analysis of composite structures, Computational Materials Science 16 (1-4) (1999) 344–354.  doi:10.1016/S0927-0256(99)00077-4.
URL https://linkinghub.elsevier.com/retrieve/pii/S0927025699000774

  * [12] K. Matouš, M. G. D. Geers, V. G. Kouznetsova, A. Gillman, A review of predictive nonlinear theories for multiscale modeling of heterogeneous materials, Journal of Computational Physics 330 (2017) 192–220.  doi:10.1016/j.jcp.2016.10.070.
URL https://linkinghub.elsevier.com/retrieve/pii/S0021999116305782

  * [13] C. Miehe, A. Koch, Computational micro-to-macro transitions of discretized microstructures undergoing small strains, Archive of Applied Mechanics (Ingenieur Archiv) 72 (4-5) (2002) 300–317.  doi:10.1007/s00419-002-0212-2.
URL http://link.springer.com/10.1007/s00419-002-0212-2

  * [14] A. Quarteroni, A. Manzoni, F. Negri, Reduced Basis Methods for Partial Differential Equations, Vol. 92 of UNITEXT, Springer International Publishing, Cham, 2016.  doi:10.1007/978-3-319-15431-2.
URL http://link.springer.com/10.1007/978-3-319-15431-2

  * [15] J. S. Hesthaven, G. Rozza, B. Stamm, Certified Reduced Basis Methods for Parametrized Partial Differential Equations, SpringerBriefs in Mathematics, Springer International Publishing, Cham, 2016.  doi:10.1007/978-3-319-22470-1.
URL http://link.springer.com/10.1007/978-3-319-22470-1

  * [16] J. Yvonnet, Q.-C. He, The reduced model multiscale method (R3M) for the non-linear homogenization of hyperelastic media at finite strains, Journal of Computational Physics 223 (1) (2007) 341–368.  doi:10.1016/j.jcp.2006.09.019.
URL https://linkinghub.elsevier.com/retrieve/pii/S0021999106004402

  * [17] A. Radermacher, B. A. Bednarcyk, B. Stier, J. Simon, L. Zhou, S. Reese, Displacement-based multiscale modeling of fiber-reinforced composites by means of proper orthogonal decomposition, Advanced Modeling and Simulation in Engineering Sciences 3 (1) (2016) 29.  doi:10.1186/s40323-016-0082-8.
URL http://amses-journal.springeropen.com/articles/10.1186/s40323-016-0082-8

  * [18] Y. Maday, N. Cuong Nguyen, A. T. Patera, S. H. Pau, A general multipurpose interpolation procedure: the magic points, Communications on Pure & Applied Analysis 8 (1) (2009) 383–404.  doi:10.3934/cpaa.2009.8.383.
URL http://aimsciences.org//article/doi/10.3934/cpaa.2009.8.383

  * [19] S. Chaturantabut, D. C. Sorensen, Nonlinear Model Reduction via Discrete Empirical Interpolation, SIAM Journal on Scientific Computing 32 (5) (2010) 2737–2764.  doi:10.1137/090766498.
URL http://epubs.siam.org/doi/10.1137/090766498

  * [20] C. Farhat, P. Avery, T. Chapman, J. Cortial, Dimensional reduction of nonlinear finite element dynamic models with finite rotations and energy-based mesh sampling and weighting for computational efficiency, International Journal for Numerical Methods in Engineering 98 (9) (2014) 625–662.  doi:10.1002/nme.4668.
URL https://onlinelibrary.wiley.com/doi/10.1002/nme.4668

  * [21] D. Ryckelynck, Hyper-reduction of mechanical models involving internal variables, International Journal for Numerical Methods in Engineering 77 (1) (2009) 75–89.  doi:10.1002/nme.2406.
URL https://onlinelibrary.wiley.com/doi/10.1002/nme.2406

  * [22] M. Yano, A. T. Patera, An LP empirical quadrature procedure for reduced basis treatment of parametrized nonlinear PDEs, Computer Methods in Applied Mechanics and Engineering 344 (2019) 1104–1123.  doi:10.1016/j.cma.2018.02.028.
URL https://linkinghub.elsevier.com/retrieve/pii/S0045782518301087

  * [23] J. Hernández, M. Caicedo, A. Ferrer, Dimensional hyper-reduction of nonlinear finite element models via empirical cubature, Computer Methods in Applied Mechanics and Engineering 313 (2017) 687–722.  doi:10.1016/j.cma.2016.10.022.
URL https://linkinghub.elsevier.com/retrieve/pii/S004578251631355X

  * [24] J. Hernández, J. Oliver, A. Huespe, M. Caicedo, J. Cante, High-performance model reduction techniques in computational multiscale homogenization, Computer Methods in Applied Mechanics and Engineering 276 (2014) 149–189.  doi:10.1016/j.cma.2014.03.011.
URL https://linkinghub.elsevier.com/retrieve/pii/S0045782514000978

  * [25] D. Soldner, B. Brands, R. Zabihyan, P. Steinmann, J. Mergheim, A numerical study of different projection-based model reduction techniques applied to computational homogenisation, Computational Mechanics 60 (4) (2017) 613–625.  doi:10.1007/s00466-017-1428-x.
URL http://link.springer.com/10.1007/s00466-017-1428-x

  * [26] A. Radermacher, S. Reese, POD-based model reduction with empirical interpolation applied to nonlinear elasticity, International Journal for Numerical Methods in Engineering 107 (6) (2016) 477–495.  doi:10.1002/nme.5177.
URL https://onlinelibrary.wiley.com/doi/10.1002/nme.5177

  * [27] M. Caicedo, J. L. Mroginski, S. Toro, M. Raschi, A. Huespe, J. Oliver, High Performance Reduced Order Modeling Techniques Based on Optimal Energy Quadrature: Application to Geometrically Non-linear Multiscale Inelastic Material Modeling, Archives of Computational Methods in Engineering 26 (4) (2019) 771–792.  doi:10.1007/s11831-018-9258-3.
URL http://link.springer.com/10.1007/s11831-018-9258-3

  * [28] M. Raschi, O. Lloberas-Valls, A. Huespe, J. Oliver, High performance reduction technique for multiscale finite element modeling (HPR-FE2): Towards industrial multiscale FE software, Computer Methods in Applied Mechanics and Engineering 375 (2021) 113580.  doi:10.1016/j.cma.2020.113580.
URL https://doi.org/10.1016/j.cma.2020.113580

  * [29] T. Guo, F. A. Silva, O. Rokoš, K. Veroy, Learning constitutive models from microstructural simulations via a non-intrusive reduced basis method: Extension to geometrical parameterizations, Computer Methods in Applied Mechanics and Engineering 401 (2022) 115636.  doi:10.1016/j.cma.2022.115636.
URL https://linkinghub.elsevier.com/retrieve/pii/S0045782522005916

  * [30] G. Rozza, D. B. P. Huynh, A. T. Patera, Reduced Basis Approximation and a Posteriori Error Estimation for Affinely Parametrized Elliptic Coercive Partial Differential Equations, Archives of Computational Methods in Engineering 15 (3) (2008) 229–275.  doi:10.1007/s11831-008-9019-9.
URL http://link.springer.com/10.1007/s11831-008-9019-9

  * [31] G. Stabile, M. Zancanaro, G. Rozza, Efficient geometrical parametrization for finite‐volume‐based reduced order methods, International Journal for Numerical Methods in Engineering 121 (12) (2020) 2655–2682.  doi:10.1002/nme.6324.
URL https://onlinelibrary.wiley.com/doi/10.1002/nme.6324

  * [32] T.-M. Yao, K. K. Choi, 3-D shape optimal design and automatic finite element regridding, International Journal for Numerical Methods in Engineering 28 (2) (1989) 369–384.  doi:10.1002/nme.1620280209.
URL https://onlinelibrary.wiley.com/doi/10.1002/nme.1620280209

  * [33] A. Cuitino, M. Ortiz, A material-independent method for extending stress update algorithms from small-strain plasticity to finite plasticity with multiplicative kinematics, Engineering computations (1992).
  * [34] J. C. Simo, T. J. R. Hughes, Computational inelasticity, Vol. 7, Springer Science & Business Media, 2006.
  * [35] M. Guo, J. S. Hesthaven, Reduced order modeling for nonlinear structural analysis using Gaussian process regression, Computer Methods in Applied Mechanics and Engineering 341 (2018) 807–826.  doi:10.1016/j.cma.2018.07.017.
URL https://linkinghub.elsevier.com/retrieve/pii/S0045782518303487

  * [36] T. Guo, O. Rokoš, K. Veroy, Learning constitutive models from microstructural simulations via a non-intrusive reduced basis method, Computer Methods in Applied Mechanics and Engineering 384 (2021) 113924.  doi:10.1016/j.cma.2021.113924.
URL
http://arxiv.org/abs/2104.04451https://linkinghub.elsevier.com/retrieve/pii/S0045782521002619http://dx.doi.org/10.1016/j.cma.2021.113924

  * [37] M. Schneider, D. Wicht, Superconvergence of the effective Cauchy stress in computational homogenization of inelastic materials, International Journal for Numerical Methods in Engineering (11 2022).  doi:10.1002/nme.7149.
URL https://onlinelibrary.wiley.com/doi/10.1002/nme.7149

  * [38] K. Bertoldi, M. C. Boyce, Mechanically triggered transformations of phononic band gaps in periodic elastomeric structures, Physical Review B 77 (5) (2 2008).  doi:10.1103/physrevb.77.052105.
  * [39] O. Rokoš, M. M. Ameen, R. H. J. Peerlings, M. G. D. Geers, Micromorphic computational homogenization for mechanical metamaterials with patterning fluctuation fields, Journal of the Mechanics and Physics of Solids 123 (2019) 119–137.  doi:10.1016/j.jmps.2018.08.019.
  * [40] Z. Zhang, A. O. Krushynska, Programmable shape-morphing of rose-shaped mechanical metamaterials, APL Materials 10 (8) (2022) 080701.  doi:10.1063/5.0099323.
URL https://aip.scitation.org/doi/10.1063/5.0099323

  * [41] M. Mozaffar, R. Bostanabad, W. Chen, K. Ehmann, J. Cao, M. A. Bessa, Deep learning predicts path-dependent plasticity, Proceedings of the National Academy of Sciences 116 (52) (2019) 26414–26420.  doi:10.1073/pnas.1911815116.
URL https://pnas.org/doi/full/10.1073/pnas.1911815116

  * [42] L. Wu, V. D. Nguyen, N. G. Kilingar, L. Noels, A recurrent neural network-accelerated multi-scale model for elasto-plastic heterogeneous materials subjected to random cyclic and non-proportional loading paths, Computer Methods in Applied Mechanics and Engineering 369 (2020) 113234.  doi:10.1016/j.cma.2020.113234.
URL https://linkinghub.elsevier.com/retrieve/pii/S0045782520304199