# Data-driven finite element computation of open-cell foam structures

Tim Fabian Korzeniowski  Kerstin Weinberg  Lehrstuhl für Festkörpermechanik,
Universität Siegen, D-57076 Siegen, Germany

###### Abstract

This paper presents a model-free data-driven strategy for linear and non-
linear finite element computations of open-cell foam. Employing sets of
material data, the data-driven problem is formulated as the minimization of a
distance function subjected to the physical constraints from the kinematics
and the balance laws of the mechanical problem.

The material data sets of the foam are deduced here from representative
microscopic material volumes. These volume elements capture the stochastic
characteristics of the open-cell microstructure and the properties of the
polyurethane material. Their computation provides the required stress-strain
data used in the macroscopic finite element computations without postulating a
specific constitutive model. The paper shows how to efficiently derive
suitable material data sets for the different (non-)linear and (an-)isotropic
material behavior cases.

Exemplarily, we compare data-driven finite element computations with
linearized and finite deformations and show that a linear kinematic is
sufficiently accurate to capture the material’s non-linearity up to 50% of
straining. The numerical example of a rubber sealing illustrates possible
areas of application, the expenditure, and the proposed strategy’s
versatility.

###### keywords:

data-driven mechanics, finite element analysis, open-cell polyurethane foam,
stochastic microstructure, representative volume elements

###### MSC:

[2020] 65M60

††journal: Journal of LaTeX Templates

##  1 Introduction

Solid foams, typically used for thermal insulation, vibration damping and
sandwich-structured composites, are lightweight engineering materials with a
heterogeneous cellular microstructure. Their design and computation require
either strong simplifications or sophisticated material models. To avoid
these, a direct use of measurement data in the computation would be desirable.
Kirchdoerfer and Ortiz [1] recently introduced such a model-free computing
method, the so-called data-driven (DD) computational mechanics, which
incorporates data directly into the finite element method (FEM), thus
bypassing the need for an explicit constitutive model. The numerical approach
was supplemented by mathematical investigations in Conti et al. [2],
demonstrating that the DD-FEM comprises the classical definition of solid
mechanics boundary value problems (BVPs). Whereas the original work [1] was
proposed in the framework of linear elasticity, several extensions to two- and
three-dimensional problems followed [3, 4, 5, 6]. Further contributions focus
on efficient data management strategies [7, 8, 9], and extend the method to
elastodynamics [4, 10], non-linear elasticity [11, 12, 13], inelasticity [14,
15], and diffusion [16]. All these contributions consider the required data to
be known.

For solid mechanic problems, the material data set which replaces the
constitutive model needs to represent strains and stresses in the dimension of
the problem. They have to encompass an almost infinite number of data tuples,
covering every possible mechanical state under consideration. Such data sets
can be gained experimentally, using in-situ computed tomography for example,
but these sorts of measurements are expensive and for many problems not
available. An alternative way of data acquisition can be the use of
computational ‘material testing’ in the sense of sample calculations for the
heterogeneous material. In the classical approach, representative volume
elements (RVEs), which allow for a detailed description of the microscopic
characteristics, are used to derive effective material properties [17]. For
the finite element analysis (FEA) of engineering components this technique
leads directly to computational homogenization and FE2 computations [18, 19],
where the two- or multi-scale problem is decoupled into nested subproblems,
i.e. a microscopic BVP and a macroscopic BVP. For the macroscopic BVP to be
solved, the required effective quantities are obtained from the solutions of
the many microscopic BVPs at the material points of the discretized continuum,
[20, 21].

For foam-like materials, homogenization strategies are chosen in [22, 23] for
example. Here, we focus specifically on an open-cell elastomeric foam. The
mechanical properties of such foams are determined by the cellular
microstructure and by the material’s properties. In our study, the foam is
made of polyurethane (PUR), a common elastomer produced by a reaction of an
isocyanate with a polyol in the presence of a catalyst and some additives.
Adding a certain amount of water to the reaction produces $\mathrm{CO_{2}}$,
which results in the formation of gaseous bubbles, [24, 25]. After
solidification, a foam with defined porosity but an irregular cellular
microstructure is ‘baked’. This variability can be captured by stochastically
generated RVEs, based on the same statistical parameters but with different
expressions of the characteristics from one realization to another.

The material data sets were generated synthetically in most previous
publications on DD-FEM [1, 11, 14]. In our contribution, we suggest utilizing
the stochastic RVEs to generate the material database required for data-driven
computations either a priori or on the fly. Exemplarily, we consider data
tuples of a deformation or strain measure $\bm{e}$ and a corresponding stress
state $\bm{s}$, which give at datum $i$ the tuple
$\bm{z}_{i}=(\bm{e},\bm{s})_{i}$. All tuples may depend on space and time; we
drop the arguments here for readability. Also, in view of finite deformation
mechanics, we leave the specific strain and stress measures open at this
point. Then, a typical material data set $\mathcal{D}$ consists of a finite
number of local states

$\displaystyle\mathcal{D}=\\{\bm{z}_{i},i=1,\dots,n\\}\,.$

The data-driven solution is then the field of states
$\mbox{\boldmath$\mathrm{z}$}\in\mathcal{D}_{g}=\mathcal{D}\times\mathcal{D}\times\mathcal{D}\times\dots$
of the global state space $\mathcal{D}_{g}$ that is closest to satisfying
kinematic compatibility and physical balance laws. Denoting the set of
compatible deformation fields in equilibrium by $\mathcal{C}$, we may
formulate the data-driven problem as

$\displaystyle\mathcal{S}=\text{argmin}\left\\{d(\mbox{\boldmath$\mathrm{z}$},\mathcal{C})\,|\,\mbox{\boldmath$\mathrm{z}$}\in\mathcal{D}_{g}\right\\}\,$
(1)

where $\mathcal{S}$ is the (approximate) solution set corresponding to the
total data set $\mathcal{D}_{g}$ and $d(\cdot,\cdot)$ is an appropriate metric
measuring the distance between the assigned data and the compatible solution.
This distance is typically measured in an energy norm.

In order to provide the data set $\mathcal{D}$, the stochastic RVE
representing a specific foam must be subjected to all sorts of deformation
states. From the solution of these microscopic BVPs, the corresponding stress
states are obtained. In such a way we generate the state space tuples
$\bm{z}_{i}$ that describes the macroscopic behavior of the heterogeneous
material.

The paper is organized as follows: In Section 2 we outline the concept of the
DD-FEM for a general, non-linear kinematic. The corresponding finite element
approximation is also described here for two different tuples of strain and
stress measures. In Section 3 we focus on the generation of the data set. We
introduce open-cell foam volumes whose microstructure is generated
stochastically. Briefly, we describe the RVE model before we discuss
strategies to generate the data set $\mathcal{D}$ for different cases of
linear or non-linear and isotropic or anisotropic material behavior. In the
following Section 4, we present three numerical examples for data-driven
computations. At first, we compare the different DD-FEM strategies for a
simple rod under tension in Section 4.1. In 4.2 and 4.3 two and three-
dimensional computations of a rubber sealing illustrate the pros and cons of
the method for engineering applications. The paper closes with a summary in
Section 5.

##  2 Governing equations

We consider a solid of domain $\Omega_{0}$ in its reference configuration and
deforming to the current configuration under the action of external body
forces $\rho_{0}\bm{B}$ and boundary tractions $\bm{T}$; the fields in
capitals refer to the reference configuration. The solid’s deformation
$\bm{\varphi}(\bm{X}):\Omega_{0}\mapsto{\rm I}\\!{\rm R}^{3}$ is completely
described by the deformation gradient

$\displaystyle\bm{F}=\nabla{\bm{\varphi}}(\bm{X})\,.$ (2)

The work conjugate stress is the first Piola-Kirchhoff tensor $\bm{P}$; the
second Piola-Kirchhoff tensor follows via $\bm{S}=\bm{F}^{-1}\bm{P}$. The
stresses fulfill the linear and the angular momentum balances

$\displaystyle\operatorname{Div}\bm{P}+\rho_{0}\bm{{B}}$
$\displaystyle=0\qquad\quad\mbox{ in }\Omega_{0}$ (3)
$\displaystyle\bm{P}^{T}\bm{F}-\bm{F}\bm{P}^{T}$ $\displaystyle=0$ (4)

and consequently is $\bm{S}=\bm{S}^{T}$. The solid is subjected to geometrical
and static boundary conditions at its boundaries $\Gamma_{0}$ and $\Gamma_{1}$
with outward unit normal $\bm{N}$,

$\displaystyle\bm{\varphi}$
$\displaystyle=\bar{\bm{\varphi}}\qquad\quad\mbox{on }\Gamma_{0}\,,$ (5)
$\displaystyle\bm{P}\,\bm{N}$ $\displaystyle={\bm{T}}\qquad\quad\mbox{on
}\Gamma_{1}\,.$ (6)

where $\Gamma_{0}\cup\Gamma_{1}=\partial\Omega_{0}$ and
$\Gamma_{0}\cap\Gamma_{1}=\emptyset$. We aim at solving this BVP in its data-
driven form (1). To do so, we first seek a proper metric of the phase space,
using hyperelastic strain-energy densities as a guide that typically have
terms of the form $|\bm{F}|^{p}$ and $|\bm{F}^{q}|$.

###  2.1 Data-driven problem in F, P

Let us start with the simplest case, namely a Neo-Hooke-like model with $p=2$
and an elastic strain-energy density as a function of the deformation
gradient,

$\displaystyle\Psi^{e}(\bm{F})=\frac{\mu}{2}\left(\bm{F}:\bm{F}-3\right)\,.$
(7)

Via a Legendre transformation, the complementary energy density as a function
of the stress tensor $\bm{P}$

$\displaystyle\Psi^{\prime}(\bm{P})$
$\displaystyle=\max_{\bm{F}\in\operatorname{\mathbb{L}\mathrm{in}}}\left[\bm{P}:\bm{F}-\frac{\mu}{2}\left(\bm{F}:\bm{F}-3\right)\right]$
$\displaystyle=\frac{1}{2\mu}\bm{P}:\bm{P}+\frac{3\mu}{2}\,.$

Thus, the distance $d(\cdot,\cdot)$ between the compatible solution and the
assigned data can be measured locally with a penalty functional of the form

$\displaystyle\bar{\Psi}(\bm{F},\bm{P})$
$\displaystyle=\min_{(\bm{F}^{\prime},\bm{P}^{\prime})\in\mathcal{D}}\left(\Psi^{e}(\bm{F}-\bm{F}^{\prime})+\Psi^{\prime}(\bm{P}-\bm{P}^{\prime})\right)$
(8)
$\displaystyle=\min_{(\bm{F}^{\prime},\bm{P}^{\prime})\in\mathcal{D}}\left(\frac{\mu_{0}}{2}(\bm{F}-\bm{F^{\prime}}):(\bm{F}-\bm{F^{\prime}})+\frac{1}{2\mu_{0}}(\bm{P}-\bm{P^{\prime}}):(\bm{P}-\bm{P^{\prime}})\right)$

where the data set
$\mathcal{D}=\\{(\bm{F},\bm{P})_{1},\dots,(\bm{F},\bm{P})_{n}\\}$ contains
tuples of values of the deformation gradient $\bm{F}$ and the first Piola-
Kirchhoff stress tensor $\bm{P}$. We presume them to be physical meaningful
and to fulfill eq. (4). The parameter $\mu$ calibrates units and magnitudes
and is purely numerical in nature; we therefore write $\mu_{0}$ in the
following. When minimizing the local penalty functional (8), the kinematics of
the deformation (2) as well as the momentum balances (3) and (4) must hold for
all admissible states $\bm{F}$ and $\bm{P}$. Equivalently, we can formulate a
global penalty functional for the whole domain or a part of it

$\displaystyle\mathcal{W}^{\text{pen}}=\int_{\Omega_{e}}\bar{\Psi}(\bm{F},\bm{P})\,\mathrm{d}\Omega\qquad\rightarrow\text{min.},$

and we may likewise demand the equilibrium conditions in global form.

We want to start with the local form (8). The data-driven BVP then reads:
$\bar{\Psi}^{*}\rightarrow\text{min.}$ with functional

$\displaystyle\bar{\Psi}^{*}=\min_{(\bm{F}^{\prime},\bm{P}^{\prime})\in\mathcal{D}}$
$\displaystyle\left(\Psi^{e}(\bm{F}-\bm{F}^{\prime})+\Psi^{\prime}(\bm{P}-\bm{P}^{\prime})\right)$
$\displaystyle+\bm{\lambda}\left(\nabla\cdot\bm{P}+\varrho_{0}\bm{B}\right)+\bm{\lambda}_{2}:\left(\bm{F}\bm{P}^{T}-\bm{P}^{T}\bm{F}\right)\,.$

The additional Lagrange multiplier field $\bm{\lambda}$ enforces the static
equilibrium (3) and is a displacement-like vector field, the tensor-valued
field $\bm{\lambda}_{2}$ enforces the angular momentum balance (4). For the
optimization of $\bar{\Psi}^{*}$ with a finite element discretization in
Lagrangian form, we formulate the deformation gradient with the local
displacements $\bm{u}$

$\displaystyle\bm{F}={\bm{I}+\nabla\bm{u}}$ (9)

so that the functional depends on $\bm{u}$, $\bm{P},\bm{\lambda}$ and
$\bm{\lambda}_{2}$. Assuming additionally that the optimal data tuples
$(\bm{F}^{*},\bm{P}^{*})$ are already determined, the functional to be
optimized reads

$\displaystyle\bar{\Psi}^{*}(\bm{u},\bm{P},\bm{\lambda},\bm{\lambda}_{2})=$
$\displaystyle\frac{\mu_{0}}{2}(\bm{I}+\nabla\bm{u}-\bm{F}^{*}):(\bm{I}+\nabla\bm{u}-\bm{F}^{*})+\frac{1}{2\mu_{0}}(\bm{P}-\bm{P}^{*}):(\bm{P}-\bm{P}^{*})$
$\displaystyle+\bm{\lambda}\cdot(\nabla\cdot{\bm{P}}+\varrho_{0}\bm{B})+\bm{\lambda}_{2}:\left((\bm{I}+\nabla\bm{u})\bm{P}^{T}-\bm{P}^{T}(\bm{I}+\nabla\bm{u})\right)$

and its total variation gives the equations:

$\displaystyle\delta_{\bm{u}}\bar{\Psi}^{*}=0:\quad\rightarrow$
$\displaystyle\mu_{0}{\nabla\delta\bm{u}}:({\bm{I}+\nabla\bm{u}}-\bm{F}^{*})=$
$\displaystyle\bm{\lambda}_{2}:(\bm{P}^{T}{\nabla\delta\bm{u}}-{\nabla\delta\bm{u}}\,\bm{P}^{T})$
(10) $\displaystyle\delta_{\bm{P}}\bar{\Psi}^{*}=0:\quad\rightarrow$
$\displaystyle\tfrac{1}{\mu_{0}}\delta\bm{P}:(\bm{P}-\bm{P}^{*})+\bm{\lambda}\cdot(\nabla\cdot\delta\bm{P})=$
$\displaystyle\bm{\lambda}_{2}:(\delta\bm{P}^{T}\,\nabla\bm{u}+\nabla\bm{u}\;\delta\bm{P}^{T})$
$\displaystyle\delta_{\bm{\lambda}}\bar{\Psi}^{*}=0:\quad\rightarrow$
$\displaystyle\delta\bm{\lambda}\cdot\left(\nabla\cdot\bm{P}+\varrho_{0}\bm{B}\right)=$
$\displaystyle 0$
$\displaystyle\delta_{\bm{\lambda}_{2}}\bar{\Psi}^{*}=0:\quad\rightarrow$
$\displaystyle\delta\bm{\lambda}_{2}:\left(\nabla\bm{u}\bm{P}^{T}-\bm{P}^{T}\nabla\bm{u}\right)=$
$\displaystyle 0$

We want to express these equations in terms of the unknown displacements
$\bm{u}$ and the the multiplier field $\bm{\lambda}$, which can conveniently
be approximated by finite elements. We start with eq. (10)4, and note that for
any $\delta\bm{\lambda}_{2}$ its total variation equals the sum of the right-
hand sides of eq. (10)1 and (10)2. Consequently, the sums of the left-hand
sides also equal. Both terms are now formulated in the integral mean over an
(element) domain $\Omega_{e}$, and the second term is re-arranged using
Gauss’s theorem and $\delta\bm{P}=0$ at the boundary. With eq. (9) we obtain

$\displaystyle\int_{\Omega_{e}}\mu_{0}{\nabla\delta\bm{u}}:(\bm{F}-\bm{F}^{*})\,\mathrm{d}\Omega=-\int_{\Omega_{e}}\frac{1}{\mu_{0}}\delta\bm{P}:(\bm{P}-\bm{P}^{*}+\mu_{0}\nabla\bm{\lambda})\,\mathrm{d}\Omega\,.$
(11)

At next we integrate eq. (10)3 over domain $\Omega_{e}$ and obtain after re-
arrangement with the boundary conditions (6)

$\displaystyle\int_{\Omega_{e}}\nabla\delta\bm{\lambda}:\bm{P}\,\mathrm{d}\Omega$
$\displaystyle=\int_{\partial\Omega_{e}}\delta\bm{\lambda}\cdot\bm{T}\,\mathrm{d}\Omega+\int_{\Omega_{e}}\delta\bm{\lambda}\cdot\varrho_{0}\bm{B}\,\mathrm{d}\Omega=:\bm{f}^{\text{ext}}\,.$
(12)

For further evaluation of eq. (12) we need an expression for the stresses
$\bm{P}$, which we want to deduce from relation (11). We remark that a
classical finite element approximation provides the displacements as the
primary solution. With the argument, that then the left-hand side of (11) is
close to zero, the bracketed term on the right-hand side of (11) must vanish
for all $\delta\bm{P}$. Thus, we derive for the stresses the expression

$\displaystyle\bm{P}=\bm{P}^{*}-\mu_{0}\nabla\bm{\lambda}$ (13)

which we now substitute into (10)3 resp. (12). What remains is

$\displaystyle\int_{\Omega_{e}}\mu_{0}\nabla\delta\bm{\lambda}:\nabla\bm{\lambda}\,\mathrm{d}\Omega=\bm{f}^{\text{ext}}-\int_{\Omega_{e}}\nabla\delta\bm{\lambda}:\bm{P}^{*}\,\mathrm{d}\Omega$
(14)

where $\bm{f}_{e}$ abbreviates the right-hand side of (12). Finally, we also
state the integral form of eq. (10)1,

$\displaystyle\int_{\Omega_{e}}\mu_{0}\nabla\delta\bm{u}:\nabla\bm{u}\,\mathrm{d}\Omega$
$\displaystyle=\int_{\Omega_{e}}\mu_{0}\nabla\delta\bm{u}:(\bm{F}^{*}-\bm{I})\,\mathrm{d}\Omega\,.$
(15)

Equations (14) and (15) formulate the data-driven problem we are looking for.
Obviously, the resulting system of equations is coupled but linear in $\bm{u}$
and $\bm{\lambda}$. The reason for this linearity is the stress expression
(13). Without this approximation it is not possible to reduce the problem to
the unknown fields $\bm{u}$ and $\bm{\lambda}$ for the energy density (7).

###  2.2 Finite element formulation in F, P

For the numerical simulation we decompose the domain $\Omega_{0}$ of dimension
$\mathrm{d}$ into finite elements $\Omega_{e}$,

$\displaystyle\Omega_{0}\approx\tilde{\Omega}_{0}=\bigcup_{\mathcal{E}}\Omega_{e}$

where ${\mathcal{E}}=\\{1,\dots,n_{el}\\}$ is the set of all elements and
$\cup$ abbreviates the assembly. We employ the usual ansatz for the unknown
fields and their variations

$\displaystyle\bm{u}(\bm{x})\approx\tilde{\bm{u}}(\bm{x})$
$\displaystyle=\mbox{\boldmath$\mathrm{N}$}\hat{\bm{u}}\qquad\qquad\delta\tilde{\bm{u}}(\bm{x})=\mbox{\boldmath$\mathrm{N}$}\,\delta\hat{\bm{u}}$
(16) $\displaystyle\bm{\lambda}(\bm{x})\approx\tilde{\bm{\lambda}}(\bm{x})$
$\displaystyle=\mbox{\boldmath$\mathrm{N}$}\hat{\bm{\lambda}}\qquad\qquad\delta\tilde{\bm{u}}(\bm{x})=\mbox{\boldmath$\mathrm{N}$}\,\delta\hat{\bm{\lambda}}$

where the matrix $\mathrm{N}$ contains the shape functions $N_{k}$ for every
degree of freedom $k=1,\dots,n_{\text{dof}}$. The vectors $\hat{\bm{u}}$ and
$\hat{\bm{\lambda}}$ denote the nodal displacements and Langrange parameters.
The gradients are summarized in matrix
$\mbox{\boldmath$\mathrm{B}$}\equiv\nabla\mbox{\boldmath$\mathrm{N}$}$; the
subscript $e$ refers to one element.

In that way we obtain the discretization of eqs. (14) and (15)

$\displaystyle\int_{\Omega_{e}}\mu_{0}{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}{\mbox{\boldmath$\mathrm{B}$}^{e}}\,\mathrm{d}\Omega\
\hat{\bm{u}}_{e}$
$\displaystyle=\int_{\Omega_{e}}\mu_{0}{\mbox{\boldmath$\mathrm{B}$}^{e}}(\mbox{\boldmath$\mathrm{F}$}^{*}-\mbox{\boldmath$\mathrm{I}$})\,\mathrm{d}\Omega$
(17)
$\displaystyle\int_{\Omega_{e}}\mu_{0}{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}{\mbox{\boldmath$\mathrm{B}$}^{e}}\,\mathrm{d}\Omega\
\hat{\bm{\lambda}}_{e}$
$\displaystyle=\mbox{\boldmath$\mathrm{f}$}^{e}-\int_{\Omega_{e}}{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}{\mbox{\boldmath$\mathrm{P}$}^{*}}\,\mathrm{d}\Omega$

with the data tuples
$(\mbox{\boldmath$\mathrm{F}$}^{*},\mbox{\boldmath$\mathrm{P}$}^{*})$ in
vector form and the elementwise force vector

$\displaystyle\mbox{\boldmath$\mathrm{f}$}^{e}=\int_{\Omega_{e}}\varrho_{0}\mbox{\boldmath$\mathrm{N}$}^{T}\mbox{\boldmath$\mathrm{B}$}\,\mathrm{d}\Omega+\int_{\Gamma_{1e}}\mbox{\boldmath$\mathrm{N}$}^{T}\mbox{\boldmath$\mathrm{T}$}\,\mathrm{d}\Gamma\,.$
(18)

The resulting finite element system is given in Appendix 1.

###  2.3 Data-driven problem in C, S

To indeed account for a non-linear kinematic of the deformation, an energy
density with a higher order term is necessary. We choose $q=4$, i.e.
$|\bm{F}^{4}|$ in density (7) and describe the deformation with the invariant
right Cauchy-Green tensor $\bm{C}=\bm{F}^{T}F$. The elastic strain-energy
density is then

$\displaystyle\Psi^{e}(\bm{C})=\frac{\mu_{0}}{2}\left(\bm{C}:\bm{C}-3\right)\,.$

with the corresponding conjugate, the second Piola-Kirchhoff stress tensor
$\bm{S}$. Again, via a Legendre transformation we obtain the complementary
energy density,

$\displaystyle\Psi^{\prime}(\bm{S})$
$\displaystyle=\frac{1}{2\mu_{0}}\bm{S}:\bm{S}+\frac{3\mu_{0}}{2}\,,$

and so the local penalty functional reads

$\displaystyle\bar{\Psi}(\bm{C},\bm{S})$
$\displaystyle=\min_{(\bm{C}^{\prime},\bm{S}^{\prime})\in\mathcal{D}}\left(\frac{\mu_{0}}{2}(\bm{C}-\bm{C^{\prime}}):(\bm{C}-\bm{C^{\prime}})+\frac{1}{2\mu_{0}}(\bm{S}-\bm{S^{\prime}}):(\bm{S}-\bm{S^{\prime}})\right)$
(19)

where the data set now contains $(\bm{C},\bm{S})$ tuples,
$\mathcal{D}=\\{(\bm{C},\bm{S})_{1},\dots,(\bm{C},\bm{S})_{n}\\}$. The
minimization of functional (19) has to be performed under the equilibrium
constraint (3), i.e. $\operatorname{Div}(\bm{F}\bm{S})+\varrho_{0}\bm{B}=0$,
the balance of angular momentum is fulfilled a priori. Expressing $\bm{C}$
with the displacement field $\bm{u}$

$\displaystyle\bm{C}={{\nabla\bm{u}}+{\nabla\bm{u}}^{T}+{\nabla\bm{u}}^{T}{\nabla\bm{u}}+\bm{I}}$

we obtain the functional to be optimized in terms of $\bm{u}$,
$\bm{S},\bm{\lambda}$. For optimal data tuples $(\bm{C}^{*},\bm{S}^{*})$ it
reads

$\displaystyle\bar{\Psi}^{*}(\bm{u},\bm{S},\bm{\lambda})=$
$\displaystyle\frac{\mu_{0}}{2}({{\nabla\bm{u}}+{\nabla\bm{u}}^{T}+{\nabla\bm{u}}^{T}{\nabla\bm{u}}+\bm{I}}-\bm{C}^{*})^{2}:\bm{I}+\frac{1}{2\mu_{0}}(\bm{S}-\bm{S}^{*})^{2}:\bm{I}$
$\displaystyle\qquad+\bm{\lambda}\cdot(\nabla\cdot{({\bm{I}+\nabla\bm{u}})\bm{S}}+\varrho_{0}\bm{B})$

and gives the variational equations:

$\displaystyle\delta_{\bm{u}}\bar{\Psi}^{*}=0:\quad\rightarrow$ $\displaystyle
2\mu_{0}\left({{\nabla\delta\bm{u}}+{\nabla\delta\bm{u}}^{T}{\nabla\bm{u}}}\right):(\nabla\bm{u}+{\nabla\bm{u}}^{T}+{\nabla\bm{u}}^{T}{\nabla\bm{u}}$
$\displaystyle+\bm{I}-\bm{C}^{*})$ (20)
$\displaystyle+\bm{\lambda}\cdot\left(\nabla\cdot{\left(\nabla\delta\bm{u}\,\bm{S}\right)}\right)$
$\displaystyle=0\qquad$
$\displaystyle\delta_{\bm{S}}\bar{\Psi}^{*}=0:\quad\rightarrow$
$\displaystyle\tfrac{1}{\mu_{0}}\delta\bm{S}:(\bm{S}-\bm{S}^{*})+\bm{\lambda}\cdot\left(\nabla\cdot{\left(({\bm{I}+\nabla\bm{u}})\delta\bm{S}\right)}\right)$
$\displaystyle=0$
$\displaystyle\delta_{\bm{\lambda}}\bar{\Psi}^{*}=0:\quad\rightarrow$
$\displaystyle\delta\bm{\lambda}\cdot\left(\nabla\cdot({\bm{I}+\nabla\bm{u}})\bm{S}+\varrho_{0}\bm{B}\right)$
$\displaystyle=0$

Again, we want to express these equations with the unknown displacements
$\bm{u}$ and the multiplier field $\bm{\lambda}$ in an integral mean over a
domain $\Omega_{e}$. From eq. (20)1 we obtain after re-arrangement

$\displaystyle
2\mu_{0}\int_{\Omega_{e}}\left({{\nabla\delta\bm{u}}+{\nabla\delta\bm{u}}^{T}{\nabla\bm{u}}}\right)$
$\displaystyle:({{\nabla\bm{u}}+{\nabla\bm{u}}^{T}+{\nabla\bm{u}}^{T}{\nabla\bm{u}}+\bm{I}}-\bm{C}^{*})\,\mathrm{d}\Omega$
(21)
$\displaystyle=\int_{\Omega_{e}}\nabla\bm{\lambda}:\nabla\delta\bm{u}\,\bm{S}\,\mathrm{d}\Omega\,.$

The eq. (20)3 becomes

$\displaystyle\int_{\Omega_{e}}\nabla\delta\bm{\lambda}\,{\left(\bm{I}+\nabla\bm{u}\right)^{T}}:\bm{S}\,\mathrm{d}\Omega=\int_{\Omega_{e}}\delta\bm{\lambda}\cdot\bm{T}\,\mathrm{d}\Omega+\int_{\Omega_{e}}\delta\bm{\lambda}\cdot\varrho_{0}\bm{B}\,\mathrm{d}\Omega$
$\displaystyle=:\bm{f}^{\text{ext}}\,,$ (22)

and from eq. (20)2 it follows

$\displaystyle\int_{\Omega_{e}}\frac{1}{\mu_{0}}(\bm{S}-\bm{S}^{*}-\mu_{0}\nabla\bm{\lambda}{\left(\bm{I}+\nabla\bm{u}\right)^{T}}):\delta\bm{S}\,\mathrm{d}\Omega$
$\displaystyle=0\,$

from where we derive the expression for the stresses,

$\displaystyle\bm{S}=\bm{S}^{*}+\mu_{0}\nabla\bm{\lambda}{\left(\bm{I}+\nabla\bm{u}\right)^{T}}$

which now depend on the displacements and the multiplier field and are clearly
non-linear. This expression is put into eq. (20)3 resp. (22), and it results

$\displaystyle\int_{\Omega_{e}}\nabla\delta\bm{\lambda}{\left(\bm{I}+\nabla\bm{u}\right)^{T}}:\left(\bm{S}^{*}+\mu_{0}\nabla\bm{\lambda}{\left(\bm{I}+\nabla\bm{u}\right)^{T}}\right)\,\mathrm{d}\Omega$
$\displaystyle=\bm{f}^{\text{ext}}$ (23)

where $\bm{f}^{\text{ext}}$ again abbreviates the right side of (22). The eqs.
(21) and (23) form the wanted coupled and non-linear system of equations for
the fields $\bm{u}$ and $\bm{\lambda}$. We summarize them again as residual
equations, writing for the sake of clarity
$\bm{F}(\bm{u})={\bm{I}+\nabla\bm{u}}$:

$\displaystyle\int_{\Omega_{e}}{\nabla\delta\bm{u}}^{T}:\bigg{(}2\mu_{0}{\bm{F}(\bm{u})^{T}}\Big{(}\bm{F}(\bm{u})^{T}\bm{F}(\bm{u})-\bm{C}^{*}\Big{)}$
(24)
$\displaystyle-\nabla\bm{\lambda}^{T}\bm{S}^{*}-\mu_{0}\nabla\bm{\lambda}^{T}\nabla\bm{\lambda}{\bm{F}(\bm{u})^{T}}\bigg{)}\,\mathrm{d}\Omega$
$\displaystyle=0$
$\displaystyle\int\limits_{\Omega_{e}}\nabla\delta\bm{\lambda}:\left(\mu_{0}{\bm{F}(\bm{u})}\nabla\bm{\lambda}{\bm{F}(\bm{u})^{T}}+{\bm{F}(\bm{u})}\bm{S}^{*}\right)\,\mathrm{d}\Omega-\bm{f}^{\text{ext}}$
$\displaystyle=0\,.$ (25)

###  2.4 Finite element formulation in C, S

The finite element discretization of the system (24-25) is performed with the
ansatz (16). The deformation gradient is calculated in each integration point
as

$\displaystyle\mbox{\boldmath$\mathrm{F}$}^{e}=\mbox{\boldmath$\mathrm{I}$}+\mbox{\boldmath$\mathrm{B}$}^{e}\hat{\bm{u}}^{e}$

where $\mbox{\boldmath$\mathrm{F}$}^{e}$ has the suitable vector form. The
discretized system of residual equations is then:

$\displaystyle{\mbox{\boldmath$\mathrm{R}$}}_{u}$
$\displaystyle=\int_{\Omega_{e}}{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}\bigg{(}2\mu_{0}{\mbox{\boldmath$\mathrm{F^{e}}$}^{T}}\Big{(}\mbox{\boldmath$\mathrm{F^{e}}$}^{T}\mbox{\boldmath$\mathrm{F^{e}}$}-\mbox{\boldmath$\mathrm{C^{*}}$}\Big{)}-{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}\hat{\bm{\lambda}}\mbox{\boldmath$\mathrm{S}$}^{*}-\mu_{0}{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}\hat{\bm{\lambda}}{\mbox{\boldmath$\mathrm{B}$}^{e}}\hat{\bm{\lambda}}{\mbox{\boldmath$\mathrm{F}$}^{e}}^{T}\bigg{)}\,\mathrm{d}\Omega=0$
(26) $\displaystyle\mbox{\boldmath$\mathrm{R}$}_{\lambda}$
$\displaystyle=\int_{\Omega_{e}}\mu_{0}{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}{\mbox{\boldmath$\mathrm{F}$}^{e}}{\mbox{\boldmath$\mathrm{B}$}^{e}}\hat{\bm{\lambda}}{\mbox{\boldmath$\mathrm{F}$}^{e}}^{T}+{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}{\mbox{\boldmath$\mathrm{F}$}^{e}}\mbox{\boldmath$\mathrm{S}$}^{*}\,\mathrm{d}\Omega-\mbox{\boldmath$\mathrm{f}$}^{e}=0$

The solution of the system (26) requires an iterative scheme, typically a
Newton-Rhapson iteration. The resulting finite element expressions are given
in Appendix 2.

##  3 Recording of the material data sets

The material data sets required for the DD-FEM are gained from systematic
computations of the representative microscopic foam volumes. Hereby we expect
the results to be physically meaningful, e.g., symmetric and regular. Also,
the specific choice of strain-stress tuples $\bm{z}_{i}$ is of minor
importance; they can be converted into each other using the common relations
of continuums mechanics. Therefore the data sets

$\displaystyle\mathcal{D}_{F}=\\{(\bm{F},\bm{P})_{i},i=1,\dots,n\\}\qquad\text{and}\qquad\mathcal{D}_{C}=\\{(\bm{C},\bm{S})_{i},i=1,\dots,n\\}\qquad$

or its small strain equivalent
$\mathcal{D}_{\epsilon}=\\{(\bm{\epsilon},\bm{\sigma})_{i},i=1,\dots,n\\}$ are
treated equally and denoted as $\mathcal{D}$ subsequently.

###  3.1 Generation of stochastic RVEs for open-cell foam

Figure 1: Procedure to generate a representative stochastic foam volume:
sphere radius and location are determined by stochastic variables. A Laguerre
tesselation assigns a domain to every sphere. Then the cells are scaled to a
given volume ratio and the skeleton is extracted afterwards.

Elastic open-cell foams are rubbery materials with a low relative density
$R=\varrho_{\text{foam}}/\varrho_{\text{elastomer}}\leq\nicefrac{{1}}{{3}}$.
Their mechanical behavior is determined by the matrix material and the
cellular microstructure, consisting of a network of ligaments connected at
junctions (vertices). An accurate description relies on the real foam’s
geometry, whose topological characteristics are gained from computed
tomography scans for example, cf. [26]. Assuming these characteristics, e.g.
the pore volume fraction, the size distribution, the coefficient of variation,
and the anisotropy factor, to be known, the following procedure is chosen to
generate the corresponding RVEs.

At first, random sphere distributions with a collective rearrangement of the
spheres are used to build a dense isotropic packing, see the left image of
Fig. 1. A force-biased packing algorithm [27] yields a very efficient
procedure to attain the desired arrangement. At next, a Laguerre tessellation,
which is a Voronoi tessellation weighted with the sphere radius [28], is used
to partition the volume into subregions, see the second image of Fig. 1. Then,
the effective density of the foam is gained by a scaling of the cells, see the
third image of Fig. 1. Finally, the edges need to be extracted. Originally
they had been shared by three cells and through the scaling they have now a
hexagonal cross section (a triangular shape with cut corners). These polygonal
struts are the foam’s ligaments and still need to get realistic cross
sections. The adjacent struts are joined in the junctions and because this
results in kinks and internal corners a smoothing algorithm is applied. We
make use of a spline-based algorithm proposed in [29]. The result can be seen
in the right image of Fig. 1.

Figure 2: RVE of a foam with 100 log-normally distributed pores, a coefficient
of variation of $0.5$, and a relative density of $0.15$ (left); its finite
element model consists of $708\,372$ tetrahedral elements (right).

For a finite element analysis the resulting RVE still needs to be meshed. We
employ here linear tetragonal solid elements; an example RVE of size
$1000\times 1000\times 1000\,$µm can be seen in Fig. 2. The displayed RVE is
meshed with 708 372 elements and a single linear computation with periodic
boundary condition takes around 10 minutes. This effort is important because
many – linear and non-linear – computations need to be conducted to obtain the
required data sets.

###  3.2 Homogenized microscopic data set

To generate a data set that describes the homogenized material behavior of the
foam, a deformation $\bar{\bm{F}}$ of the RVE is prescribed and the stress
field is computed with a linear or non-linear FEM. The homogenized stresses
$\bar{\bm{P}}$ or $\bar{\bm{S}}$ are derived by

$\displaystyle\bar{\bm{S}}$
$\displaystyle=\frac{1}{|V|}\int_{V}\bm{S}(\bm{X})\,\mathrm{d}V\,,$

i.e. for all integration points the stress components multiplied with the
corresponding volume are summarized and divided by the RVE’s volume $V$. The
data sets, which are used as an input for the DD-FEM, are then the collections
of all $(\bar{\bm{F}},\bar{\bm{P}})$ tuples or, with
$\bar{\bm{C}}=\bar{\bm{F}}^{T}\bar{\bm{F}}$ all $(\bar{\bm{C}},\bar{\bm{S}})$
tuples. Naturally, the question arises of how many computations need to be
conducted to obtain a data set that sufficiently describes a foam’s material
behavior. Our experience suggests a minimum number of $n=50-100$ data tuples
per stress component, which results in the computation of $n^{6}$ tuples for
three dimensions at least. In the following, we discuss how we gain these data
collections efficiently in cases of different material characteristics.

####  Case A: non-linear and anisotropic material

For the general situation of a direction dependent and non-linear material
response, e.g. a viscoelastic foam with elongated pores in one direction,
simplifications are difficult. Here only a brute force approach of one RVE
computation for every $(\bar{\bm{F}},\bar{\bm{P}})$ or
$(\bar{\bm{C}},\bar{\bm{S}})$ tuple gives the wanted data set $\mathcal{D}$.
Such a strategy is hardly feasible with the RVEs described above, only a
simpler beam discretization may help. Additionally, it needs to be considered
whether these volume elements are still representative.

####  Case B: non-linear and isotropic material

For an isotropic material, i.e. when the microscopic response is non-linear
but does not depend on the loading direction, we gain the possibility to
rotate the RVE’s coordinate system, see Fig. 3. Our strategy is here to sample
discrete values of the principal stretches $\lambda_{\alpha}$,
$\alpha=1\dots$d, and then transform the resulting principal stresses,

$\displaystyle\bar{\bm{S}}=\bm{Q}\bar{\bm{S}}^{\lambda_{\alpha}}\bm{Q}^{T}$

where $\bm{Q}\in\operatorname{SO(3)}3$ is the corresponding three-dimensional
rotation tensor composed of the rotations around the coordinate axes, [1].
With a stepwise rotation around these axes all deformation states can be
mapped and summarized in $\mathcal{D}$. For $n^{6}$ data then $n^{3}$
computations are needed.

Figure 3: Elongation in principle stretches and rotation to the reference
coordinate system for d=2

####  Case C: linear and anisotropic material

The computation eases significantly when the applied deformation is small and
the microscopic response is linear. Without isotropy of the material still RVE
computations with loads in every distinct direction need to be performed but
the results now can be superposed

$\displaystyle(\bar{\bm{C}},\bar{\bm{S}})_{k}=a(\bar{\bm{C}},\bar{\bm{S}})_{i}+b(\bar{\bm{C}},\bar{\bm{S}})_{j}$

with $a,b\in{\rm I}\\!{\rm R}$. Conveniently, unit loads as displayed in Fig.
4 are computed and evaluated. This strategy leads to a enormous decrease of
computational effort and will later be employed for the small deformation
data, see A.3.

Figure 4: Unit load deformations of an RVE for direction dependent linear
material.

####  Case D: linear and isotropic material

In the simplest case of isotropic, linear (elastic) material the loading
scenarios simplify even more. The material’s isotropy makes it possible to
describe the material by only two states, namely one of the three elongations
and one of the three shear states from Fig. 4. All other data tuples can be
gained by linear combinations and rotations of the coordinate system.

Here we remark that for these simplifications we assume the RVEs to give the
same averaged response, i.e., we neglect the noise induced by the stochastic
generation. The uncertainty in the data will be subject of a subsequent work.

###  3.3 Multi-level method

The computational effort of a data-driven multiscale analysis is enormous. As
outlined above, for a non-linear microscopic material response many RVEs need
to be computed to gain the required data. Additionally, data-driven
computations are numerically costly due to the search for the optimal data
points and, furthermore, for a macroscopic non-linear kinematic the solution
requires an iterative solution procedure. To reduce these computational costs,
we apply here a multi-level method which was introduced in [7]. In this
procedure we compute the solution on smaller subsets of the original set
$\mathcal{D}_{l}\subset\mathcal{D}$. After a computation at level $l$ with
data set $\mathcal{D}_{l}$, relevant data from the total set $\mathcal{D}$ are
added to the next set $\mathcal{D}_{l+1}$. Here we go a step further than in
[7], and use the multi-level method to identify regions where more RVE
calculations are needed and then calculate additional deformation states on
the fly.

Figure 5: Illustration of the multi-level method for a simple state space: The
black data points symbolize the input set of the initial level
$\mathcal{D}_{0}$; the crosses are the data points actually assigned to the
material points. The green and red points mark the data of the first and
second refinement level, $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$,
respectively.

Fig. 5 shows exemplary three data sets on a (simplified) state space. The aim
is to cover the space of possible states (black square) as good as possible by
data tuples (dots). Instead of a high density of data, at first only a coarse
set $\mathcal{D}_{0}$ is created through RVE calculations (black dots, level
0). This set $\mathcal{D}_{0}$ is used for a first DD-FEM of the structure. In
this computation, the approximate data of relevance are identified by the data
tuples which are assigned to the material points (black $\times$). Then the
data set is adaptively refined and the data-driven FEM is repeated with a set
$\mathcal{D}_{1}$. In Fig. 5 this is exemplarily illustrated by the two green
$\times$, marking the selected data of set $\mathcal{D}_{1}$. Around those two
assigned data tuples the data grid is refined. This leads to the data set
$\mathcal{D}_{2}$ (red dots). This procedure can be repeated as often as
needed to achieve a sufficient accuracy. The adaptive computational strategy
for the multi-scale simulation is illustrated in Fig. 6.

Figure 6: Illustration of the multi-level approach. A coarse data set is
generated by RVE simulations first. A macroscopic data-driven simulation
afterwards identifies region where more precise data is needed. Additional RVE
simulations are conducted to generate new data in the desired regions. The
latter two steps can be repeated as often as desired.

##  4 Numerical Examples

To demonstrate the different DD-FEM approaches, we start our computations with
a parametric study of a simple rod. Then we show the applicability of the
method for the computation of engineering components by means of a
polyurethane rubber sealing.

###  4.1 Example 1: Rod under tension

At first we investigate the influence of different data sets together with a
linear or a non-linear kinematic. For these parametric studies we consider a
one-dimensional rod of length $l=100\,$mm and cross-section $A=1\,$mm2 under
uniaxial tension. Assuming incompressibility of the material, the state of
deformation is given by

$\displaystyle\bm{F}=\begin{pmatrix}\lambda_{1}&0&0\\\
0&\lambda_{1}^{-1/2}&0\\\ 0&0&\lambda_{1}^{-1/2}\end{pmatrix},$

which means we can uniquely describe it with the lateral stretch
$\lambda_{1}$. In this example, we do not obtain the material data sets from
any particular foam RVE because this would not be incompressible. Instead we
generate the data artificially. These data sets are then used in the two
algorithms of Section 2.

#### Artificial material data generation

We consider the lateral stretch $\lambda_{1}$ to be given and assign the
corresponding first Piola-Kirchhoff stress $P_{1}$. Three sets of data are
generated, assuming one proportional stretch-stress relation and two different
strain energy densities, whereby in the latter cases we assume the continuum
mechanics relation $\bm{S}=\bm{F}^{-1}\bm{P}$ to be valid.

  * (i)

Linear data: Assuming proportionality between $\lambda_{1}$ and $P_{1}$ as
well as $\lambda_{1}^{2}$ and $S_{1}$ we gain the following data sets:

$\displaystyle(F,P)_{i}$
$\displaystyle=\left(\lambda_{1},\,c_{1}\lambda_{1}\right)_{i}$
$\displaystyle(C,S)_{i}$
$\displaystyle=\left(\lambda_{1}^{2},\,c_{1}\lambda_{1}^{2}\right)_{i}\qquad\text{for
}i=1,\dots,n$

with material constant $c_{1}=1\,$MPa.

  * (ii)

Neo-Hookean-like data: Let the data derive from an incompressible material’s
strain energy density of the form

$\displaystyle W$
$\displaystyle=c_{1}(\operatorname{tr}(\bm{F}^{T}\bm{F})-3)+p(\det\bm{F}-1)$

where $p$ is a Lagrange multiplier enforcing incompressibility, we deduce the
two data sets from the corresponding analytical solution, cf. [30].

$\displaystyle(F,P)_{i}$
$\displaystyle=\left(\lambda_{1},\,2c_{1}(\lambda_{1}-\lambda_{1}^{-2})\right)_{i}$
$\displaystyle(C,S)_{i}$
$\displaystyle=\left(\lambda_{1}^{2},\,2c_{1}(1-\lambda_{1}^{-3})\right)_{i}\qquad\text{for
}i=1,\dots,n$

The single material constant is here $c_{1}=\nicefrac{{1}}{{6}}\,$MPa.

  * (iii)

Yeoh-like data: Using a strongly non-linear strain energy density of the form

$\displaystyle W$
$\displaystyle=c_{1}(\operatorname{tr}(\bm{F}^{T}\bm{F})-3)+c_{3}(\operatorname{tr}(\bm{F}^{T}\bm{F})-3)^{3}+p(\det\bm{F}-1)$

we derive in the same way the data tuples, cf. [31]

$\displaystyle(F,P)_{i}$
$\displaystyle=\left(\lambda_{1},\,2(\lambda_{1}-\lambda_{1}^{-2})\left(c_{1}+3c_{3}(\lambda_{1}^{2}+2\lambda_{1}^{-1}-3)^{2}\right)\right)_{i}$
$\displaystyle(C,S)_{i}$
$\displaystyle=\left(\lambda_{1}^{2},\,2(1-\lambda_{1}^{-3})\left(c_{1}+3c_{3}(\lambda_{1}^{2}+2\lambda_{1}^{-1}-3)^{2}\right)\right)_{i}\quad\text{for
}i=1,\dots,n$

with constants $c_{1}=\nicefrac{{1}}{{6}}\,$MPa and
$c_{3}=\nicefrac{{1}}{{1000}}\,$MPa.

The data of the three sets are plotted in Fig. 7. Because we chose very dense
data sets with $n=10\,000$, single data points are not recognizable in the
plot but the non-linearity can clearly be seen.

Figure 7: Data of the three sets for the one-dimensional study.

#### Comparative study

We are interested in the force-displacement relation for a homogenously
elongated bar and evaluate the maximum displacement at the bars end
u$=(\lambda_{1}-1)l$. In Fig. 8 the corresponding force-displacement curves
are plotted for the proportional material data of set (i). Clearly, computing
the DD-FEM with the linear kinematic of Section 2.1 we obtain a linear u-f-
relation, whereas with the non-linear kinematic of Section 2.3 the maximum
displacement is reduced. The same situation is displayed in Fig. 9 for the DD-
FEM computation with the sub-linear Neo-Hookean set (ii) and in Fig. 10 for
the strongly non-linear data set (iii).

Figure 8: Load-displacement curves for proportional material sets (i) computed
with the linear $(\bm{P},\bm{F})$ formulation and the non-linear
$(\bm{C},\bm{S})$ formulation.

Figure 9: Load-displacement curves for the non-linear material set (ii)
computed with the linear $(\bm{P},\bm{F})$ formulation and the non-linear
$(\bm{C},\bm{S})$ formulation; two different regimes with $500\%$ straining
(left) and $100\%$ straining (right) are shown.

Figure 10: Load-displacement curves for the strongly non-linear material set
(iii) computed with the linear $(\bm{P},\bm{F})$ formulation and the non-
linear $(\bm{C},\bm{S})$ formulation; two different regimes with $100\%$
straining (left) and $60\%$ straining (right) are shown.

In case (ii) two effects work against each other: The non-linear data are
softer for larger displacements but the non-linear kinematics stiffens the
model. In other words, while the kinematic gives a higher stretch at the same
displacement, the proportionality of model (ii) gives a lower stress at the
same stretch. The difference between the two approaches is about $14\%$,
$21\%$ and $41\%$ at $\lambda_{1}=1.5$, $\lambda_{1}=2$ and $\lambda_{1}=3$,
respectively.

In Fig. 10 the strong non-linearity of the material dominates the situation.
The u-f-curves are similar but the two effects almost balance. The linear or
non-linear kinematic hardly makes a difference here because the displacements
are also significantly smaller than before. The relative difference between
both DD-FEM computations is $8\%$ and $19\%$ at $\lambda_{1}=1.5$ and
$\lambda_{1}=2$, respectively.

Summarizing we state that in the DD-FEM the structure of the material data set
dominates the solution. For moderate straining and considering some expected
noise in the available data sets, a non-linear material behavior can be mapped
with a linear kinematic sufficiently accurate. This result is important
because a non-linear kinematic with a Newton-Raphson iteration and a data
search in every step as described in A.2 is by several factors more expensive
than a linear DD-FEM computation.

###  4.2 Example 2: A rubber sealing in a plane strain state

Now we consider a typical engineering application, namely a car door sealing
made of a foamy rubber material. The computed geometry is displayed in Fig.
11.

Figure 11: Geometry of the rubber sealing

We start with a somewhat artificial situation of a plane approximation and
two-dimensional foam data. This corresponding RVE is a simple $7\times 7\,$mm2
square with 4 regularly arranged pores of radius $2\,$mm and meshed with 2601
finite elements. The RVE matrix is a Neo-Hookean material extended to the
compressible range,
$W=c_{1}(\operatorname{tr}(\bm{C})-3)+\nicefrac{{1}}{{D_{1}}}(J-1)^{2}$, with
parameters of polyurethane, $c_{1}=3.85\,$MPa and $D_{1}=0.12$ MPa-1. The
computation is performed with the commercial program Abaqus.

#### Macroscopic finite element model

The mesh of the sealing, its loading and boundary conditions are shown on
left-hand side of Fig. 12. The component is fixed in the $z$-direction. For
the distributed surface load $p$ four different values are used, namely
$p=3\,$kPa, $p=10\,$kPa, $p=20\,$kPa and $p=30\,$kPa. For comparison, linear
(material and kinematic) and non-linear (material and kinematic) finite
element analyses based on the classical material models were also performed.

#### Generation of the material data set

For data generation the principal stretches $\lambda_{1}$ and $\lambda_{2}$
were varied in 51 steps between $0.8\dots 1.2$ which lead to 2601 non-linear
finite element computations of the RVE. The resulting data set has $2601\cdot
73=189\,873$ data tuples which provide in our data-driven simulation the
initial set $\mathcal{D}_{0}$ of the multi-level approach.

Figure 12: Mesh and boundary conditions of the rubber sealing on the left
side. On the right side the load-displacement curve is displayed for the three
simulations.

#### Data-driven finite element analysis

The data-driven computations were run on three levels of refined data sets. In
Table 1 the size of the input data set and the number of data tuples which are
actually assigned to a material point (solution set $\mathcal{S}$) are
displayed for the simulation with the highest load $p=30\,$kPa.

On the right-hand side of Fig. 12 the load-displacement curves are shown for
the data-driven solution and the two classic computations. Recorded is the
displacement of the center top nodes where we have the largest displacement.
The data-driven solution, computed with the non-linear material data but the
linearized $(\bm{F},\bm{P})$ formulation of Section 2.1, gives a somewhat
stiffer response than the fully non-linear Neo-Hookean model. The difference
in the maximal displacement is about $15\%$ for moderate straining. As
outlined in the previous example, a fully non-linear kinematic of the DD-FEM
comes at the price of a much higher computational effort and so we consider
the linearized data-driven response here to describe the deformation of the
structure sufficiently.

In Fig. 13 the horizontal Cauchy stress component $\sigma_{x}$ is displayed
for $p=20\,$kPa and a good agreement of the classic FEM solution and the data-
driven solution (level 2) can be seen.

Level $l$ | data tuples in $\mathcal{D}_{l}$ | data tuples in $\mathcal{S}_{l}$
---|---|---
0 | 189.873 | 2.972
1 | 80.244 | 3.197
2 | 86.319 | 5.201
Table 1: Sizes of the input data and the number of data tuples describing the
material of the multi level computation.  Figure 13: Horizontal Cauchy stress
$\sigma_{x}$ for a classic linear FEM solution (top) and the DD-FEM (bottom)
for a surface load of $p=20\,$kPa.

###  4.3 Example 3: Three-dimensional computation of the rubber sealing

We proceed with the rubber sealing component of Section 4.2 but compute now a
fully three-dimensional deformation state with a different loading scenario.
The material is described with microscopic RVEs of polyurethane foam.

#### Macroscopic finite element model

The geometry and boundary conditions are the same as before, see Fig. 11, but
the component has now a thickness of $5\,$mm. Additional to the surface load
$p_{1}$ from atop we add two surface loads $p_{2}$ pointing into the
transversal $z$-direction and thus twisting the component, see Fig. 14. We set
$p_{1}=-7.5\,$kPa and $p_{2}=\pm 25\,$kPa.

Figure 14: Mesh and boundary conditions of the rubber sealing in Section 4.3.
An out-of plane surface load in $z$-direction is applied (red); the elements
effected by the load $p_{1}$ are in darker green.

#### Generation of the material data set

To mimic the polyurethane foam we use the RVE presented in Section 3.1 and
compute it with the material data of Section 4.3. Under the loads of Fig. 14
we expect only moderate deformations and because of that we consider a linear
behavior. Due to the stochastic random sphere packing of the RVE, small
deviations in the directions may still occur. Therefore, we chose the strategy
outlined in Case C and conduct FEAs for the six unit loads of Fig. 4. The six
evaluated data tuples are listed in the appendix. We observe only small
deviations of around $5\%$ in the different directions; they are forwarded
into the data set. Fig. 15 shows the RVE once sheared and once elongated in
$x$ direction by $\epsilon=2\%$ and the relative low stresses. The material
data set comprises stresses between $-0.2\,$MPa and $0.2\,$MPa. The number of
data tuples varies due to the multi-level approach and is listed for the
different levels of refinement in Table 2.

Figure 15: Two of the six loading scenarios of case C: On the left side a shear deformation is displayed and on the right side the RVE that is introduced in section 3.1 is elongated along the x-axis. Level $l$ | data tuples in $\mathcal{D}_{l}$ | data tuples in $\mathcal{S}_{l}$
---|---|---
0 | 1.771.561 | 2.179
1 | 1.588.491 | 8.889
2 | 6.480.081 | 12.541
3 | 9.142.389 | 14.152
4 | 10.316.808 | 14.177
Table 2: Size of the material data set $\mathcal{D}$ and number of data tuples
describing the solution $\mathcal{S}$ in the multi level computation.

#### Data-driven finite element analysis

For the data-driven computation we used the multi-level method with 4 levels
of refinement. The data refinement is stopped then because the increase of
additional data points assigned to a material point is small in the last step.
In Fig. 16 the horizontal stress $\sigma_{x}$ is plotted for the DD-FEM. The
displayed stress distribution shows that the results are satisfying for
$\sigma_{x}$.

Figure 16: View of the normal stress $\sigma_{x}$ from above in the rubber
sealing of the DD-FEM (bottom).

The corresponding shear stress $\tau_{xz}$ is displayed in Fig. 17 in the
deformed sealing. The displacement is magnified here by a factor of five so
that the deformation of the component can nicely be seen. The numerical
generation of the data can be useful for smaller stress components. In [6] it
is noted that for an equidistant data mesh smaller stress components could not
be displayed adequately. The material density was derived by the maximal
stress component and the number of data points. Since the data is generated
numerically here, one would not be bound to a minimum recording resolution of
an experimental setup and could still refine individual components here, which
was not overused however. Finally, the data-driven method manages to perform
the macroscopic FEM simulation satisfactorily.

Figure 17: Shear stress $\tau_{xz}$ in the rubber sealing for the DD-FEM
(level 4). The numerical generation allows to generate suitable data in all
directions.

##  5 Conclusion

In this paper, we propose a strategy to provide the necessary material
database for linear and non-linear data-driven finite element computations.
Typically, such data are presumed to be available from experimental
investigations but here we suggest a computational material testing instead.
With representative material volumes the micro-scale problem is investigated
and the derived homogenized data give the input for the solution of the macro-
scale problem.

Specifically we consider the DD-FEM of components made of open-cell foam. The
representative volume elements map the foam’s characteristics such as pore
volume, pore distribution, pore and ligament geometry which are usually gained
from CT scans. Their computation, in conjunction with simplifications made
possible by the spatial material behavior, enable us to generate the material
data input set.

We showed the capabilities of the methodology with three different examples.
The first example is a study of how the data-driven method behaves with non-
linear data and a linear setting and the other way round. In the two rubber
sealing examples we focus onto the material generation and simulation of the
component. From our point of view, the DD-FEM and the data-based generation of
the RVE can really excel the computation of such problems. With the methods
presented here, it is possible to numerically calculate the properties of the
foam, and the response of a component made of it, only with the help of CT
scans. No further practical or experimental action is necessary. In
particular, the DD method can help in the future design of new types of foam.
If a new foam is to be designed, its effect on the final product can be
determined in advance by varying individual target parameters in the creation
of the RVE. By selecting the target parameters in the RVE, e.g. pore volume or
pore size, it is possible to determine which microscopic properties the foam
has and how a component made of it would react to a given load. This whole
process would be purely numerical and can save costs as well as time.

Furthermore, the data-driven method can be used to insert material properties
that were not inserted from a theoretical point of view, for example non-
linear data in a linear simulation or anisotropy in an isotropic calculation.
From our point of view, data-driven FEM is an interesting and developable
methodology that should be given further attention in the future.

## Acknowledgement

We gratefully acknowledge the support of the Deutsche Forschungsgemeinschaft
(DFG) under the project WE2525/10-2 as part of the Priority Program SPP1886
“Polymorphic uncertainty modelling for the numerical design of structures”
(no. 273721697).

## Declaration of interest

The authors declare that there is no conflict of interest regarding the
publication of this article.

## References

  * Kirchdoerfer and Ortiz [2016] T. Kirchdoerfer, M. Ortiz,  Data-driven computational mechanics,  Computer Methods in Applied Mechanics and Engineering 304 (2016) 81–101.
  * Conti et al. [2018] S. Conti, S. Müller, M. Ortiz,  Data-driven problems in elasticity,  Archive for Rational Mechanics and Analysis (2018) 1–45.
  * Kirchdoerfer and Ortiz [2017] T. Kirchdoerfer, M. Ortiz,  Data driven computing with noisy material data sets,  Computer Methods in Applied Mechanics and Engineering 326 (2017) 622–641.
  * Kirchdoerfer and Ortiz [2018] T. Kirchdoerfer, M. Ortiz,  Data-driven computing in dynamics,  International Journal for Numerical Methods in Engineering 113 (2018) 1697–1710.
  * Leygue et al. [2017] A. Leygue, M. Coret, J. Réthoré, L. Stainier, E. Verron,  Data driven constitutive identification,  preprint hal-01452494v2 (2017).
  * Korzeniowski and Weinberg [2019] T. F. Korzeniowski, K. Weinberg,  A comparison of stochastic and data-driven fem approaches to problems with insufficient material data,  Computer Methods in Applied Mechanics and Engineering 350 (2019) 554–570.
  * Korzeniowski and Weinberg [2021] T. F. Korzeniowski, K. Weinberg,  A multi-level method for data-driven finite element computations,  Computer Methods in Applied Mechanics and Engineering 379 (2021) 113740\.
  * Eggersmann et al. [2021a] R. Eggersmann, L. Stainier, M. Ortiz, S. Reese,  Model-free data-driven computational mechanics enhanced by tensor voting,  Computer Methods in Applied Mechanics and Engineering 373 (2021a) 113499\.
  * Eggersmann et al. [2021b] R. Eggersmann, L. Stainier, M. Ortiz, S. Reese,  Efficient data structures for model-free data-driven computational mechanics,  Computer Methods in Applied Mechanics and Engineering 382 (2021b) 113855\.
  * Nguyen et al. [2020] L. T. K. Nguyen, M. Rambausek, M.-A. Keip,  Variational framework for distance-minimizing method in data-driven computational mechanics,  Computer Methods in Applied Mechanics and Engineering 365 (2020) 112898\.
  * Nguyen and Keip [2018] L. T. K. Nguyen, M.-A. Keip,  A data-driven approach to nonlinear elasticity,  Computers & Structures 194 (2018) 97–115.
  * Ibanez et al. [2018] R. Ibanez, E. Abisset-Chavanne, J. V. Aguado, D. Gonzalez, E. Cueto, F. Chinesta,  A manifold learning approach to data-driven computational elasticity and inelasticity,  Archives of Computational Methods in Engineering 25 (2018) 47–57.
  * Platzer [2020] A. Platzer, Finite strain data-driven computational mechanics.: From tailored data to adaptive solvers for multiscale simulations, Ph.D. thesis, École centrale de Nantes, 2020\.
  * Eggersmann et al. [2019] R. Eggersmann, T. Kirchdoerfer, S. Reese, L. Stainier, M. Ortiz,  Model-free data-driven inelasticity,  Computer Methods in Applied Mechanics and Engineering 350 (2019) 81–99.
  * Ibáñez et al. [2019] R. Ibáñez, E. Abisset-Chavanne, D. González, J.-L. Duval, E. Cueto, F. Chinesta,  Hybrid constitutive modeling: data-driven learning of corrections to plasticity models,  International Journal of Material Forming 12 (2019) 717–725.
  * Korzeniowski and Weinberg [2021] T. F. Korzeniowski, K. Weinberg,  Data-driven diffusion with uncertainties,  PAMM 20 (2021) e202000325.
  * Allaire [1992] G. Allaire,  Homogenization and two-scale convergence,  SIAM Journal on Mathematical Analysis 23 (1992) 1482–1518.
  * Miehe et al. [1999] C. Miehe, J. Schröder, J. Schotte,  Computational homogenization analysis in finite plasticity simulation of texture development in polycrystalline materials,  Computer methods in applied mechanics and engineering 171 (1999) 387–418.
  * Feyel and Chaboche [2000] F. Feyel, J.-L. Chaboche,  Fe2 multiscale approach for modelling the elastoviscoplastic behaviour of long fibre sic/ti composite materials,  Computer methods in applied mechanics and engineering 183 (2000) 309–330.
  * Abdulle et al. [2012] A. Abdulle, E. Weinan, B. Engquist, E. Vanden-Eijnden,  The heterogeneous multiscale method,  Acta Numerica 21 (2012) 1–87.
  * Cioranescu et al. [2008] D. Cioranescu, A. Damlamian, G. Griso,  The periodic unfolding method in homogenization,  SIAM Journal on Mathematical Analysis 40 (2008) 1585–1620.
  * Kanaun and Tkachenko [2007] S. Kanaun, O. Tkachenko,  Representative volume element and effective elastic properties of open cell foam materials with random microstructures,  Journal of Mechanics of Materials and Structures 2 (2007) 1607–1628.
  * Marvi-Mashhadi et al. [2018] M. Marvi-Mashhadi, C. Lopes, J. LLorca,  Modelling of the mechanical behavior of polyurethane foams by means of micromechanical characterization and computational homogenization,  International Journal of Solids and Structures 146 (2018) 154–166.
  * Kaiser [2015] W. Kaiser, Kunststoffchemie für Ingenieure: von der Synthese bis zur Anwendung, Carl Hanser Verlag GmbH Co KG, 2015.
  * Bayer [1947] O. Bayer,  Das di-isocyanat-polyadditionsverfahren (polyurethane),  Angewandte Chemie 59 (1947) 257–272.
  * Buchen et al. [2022] S. Buchen, L. Bogunia, S. Siebert, K. Weinberg,  Microstructure characterization of foamed rubber based on intrinsic volumes of $\mu$ct-images,  in preperation (2022).
  * Bargieł and Mościński [1991] M. Bargieł, J. Mościński,  C-language program for the irregular close packing of hard spheres,  Computer Physics Communications 64 (1991) 183–192.
  * Liebscher [2015] A. Liebscher,  Laguerre approximation of random foams,  Philosophical Magazine 95 (2015) 2777–2792.
  * Loop [1994] C. Loop,  Smooth spline surfaces over irregular meshes,  in: Proceedings of the 21st annual conference on Computer graphics and interactive techniques, 1994, pp. 303–310.
  * Reppel et al. [2013] T. Reppel, T. Dally, K. Weinberg,  On the elastic modeling of highly extensible polyurea,  Technische Mechanik-European Journal of Engineering Mechanics 33 (2013) 19–33.
  * Yeoh [1990] O. H. Yeoh,  Characterization of elastic properties of carbon-black-filled rubber vulcanizates,  Rubber chemistry and technology 63 (1990) 792–805.

##  A Appendix

###  A.1 Systems of equations for the F,P formulation

The linear finite element system of Eq. (17) and (18) reads:

$\displaystyle\mbox{\boldmath$\mathrm{K}$}^{e}_{u}{\hat{\bm{u}}_{e}}=\mbox{\boldmath$\mathrm{f}$}^{e}_{u}\
:\qquad$
$\displaystyle\mbox{\boldmath$\mathrm{K}$}_{u}^{e}=\mu_{0}\int_{\Omega^{e}}{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}\mbox{\boldmath$\mathrm{B}$}^{e}\,\mathrm{d}\Omega$
$\displaystyle\mbox{\boldmath$\mathrm{K}$}_{u}=\bigcup_{\mathcal{E}}\mbox{\boldmath$\mathrm{K}$}^{e}\qquad\qquad$
$\displaystyle\mbox{\boldmath$\mathrm{f}$}^{e}_{u}=\mu_{0}\int_{\Omega^{e}}{\mbox{\boldmath$\mathrm{B}$}^{e}}{\bm{\epsilon}^{*}}\,\mathrm{d}\Omega$
$\displaystyle\mbox{\boldmath$\mathrm{f}$}_{u}=\bigcup_{\mathcal{E}}\mbox{\boldmath$\mathrm{f}$}^{e}\qquad\qquad\
$
$\displaystyle\mbox{\boldmath$\mathrm{K}$}^{e}_{\lambda}\hat{\bm{\lambda}}_{e}=\mbox{\boldmath$\mathrm{f}$}^{e}_{\lambda}\
:\qquad$
$\displaystyle\mbox{\boldmath$\mathrm{K}$}_{\lambda}^{e}=\mu_{0}\int_{\Omega^{e}}{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}\mbox{\boldmath$\mathrm{B}$}^{e}\,\mathrm{d}\Omega$
$\displaystyle\mbox{\boldmath$\mathrm{K}$}_{\lambda}=\bigcup_{\mathcal{E}}\mbox{\boldmath$\mathrm{K}$}^{e}\qquad\qquad$
$\displaystyle\mbox{\boldmath$\mathrm{f}$}^{e}_{\lambda}=\mbox{\boldmath$\mathrm{f}$}^{e}-\int_{\Omega^{e}}{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}{\bm{\sigma}^{*}}\,\mathrm{d}\Omega$
$\displaystyle\mbox{\boldmath$\mathrm{f}$}_{\lambda}=\bigcup_{\mathcal{E}}\mbox{\boldmath$\mathrm{f}$}^{e}\qquad\qquad\
$

###  A.2 Systems of equations for the C,S formulation

The non-linear finite element system of Eq. 26 reads as residual equations

$\displaystyle{\mbox{\boldmath$\mathrm{R}$}}_{u}$
$\displaystyle=\int_{\Omega_{e}}{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}\bigg{(}2\mu_{0}{\mbox{\boldmath$\mathrm{F}$}^{e}}^{T}\Big{(}{\mbox{\boldmath$\mathrm{F}$}^{e}}^{T}\mbox{\boldmath$\mathrm{F}$}^{e}-\mbox{\boldmath$\mathrm{C^{*}}$}\Big{)}-{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}\hat{\bm{\lambda}}\mbox{\boldmath$\mathrm{S}$}^{*}-\mu_{0}{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}\hat{\bm{\lambda}}{\mbox{\boldmath$\mathrm{B}$}^{e}}\hat{\bm{\lambda}}{\mbox{\boldmath$\mathrm{F}$}^{e}}^{T}\bigg{)}\,\mathrm{d}\Omega=0$
$\displaystyle\mbox{\boldmath$\mathrm{R}$}_{\lambda}$
$\displaystyle=\int_{\Omega_{e}}\mu_{0}{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}{\mbox{\boldmath$\mathrm{F}$}^{e}}{\mbox{\boldmath$\mathrm{B}$}^{e}}\hat{\bm{\lambda}}{\mbox{\boldmath$\mathrm{F}$}^{e}}^{T}+{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}{\mbox{\boldmath$\mathrm{F}$}^{e}}\mbox{\boldmath$\mathrm{S}$}^{*}\,\mathrm{d}\Omega-\mbox{\boldmath$\mathrm{f}$}^{\text{ext}}=0$

A corresponding Newton-Raphson iteration step has the form

$\displaystyle\begin{bmatrix}\hat{\bm{u}}_{j+1}\\\
\hat{\bm{\lambda}}_{j+1}\end{bmatrix}=\begin{bmatrix}\hat{\bm{u}}_{j}\\\
\hat{\bm{\lambda}}_{j}\end{bmatrix}+\begin{bmatrix}\Delta\hat{\bm{u}}\\\
\Delta\hat{\bm{\lambda}}\end{bmatrix}\quad\text{ mit
}\quad\begin{bmatrix}\mbox{\boldmath$\mathrm{R}$}_{u}(\hat{\bm{u}}_{j},\hat{\bm{\lambda}}_{j})\\\
\mbox{\boldmath$\mathrm{R}$}_{\lambda}(\hat{\bm{u}}_{j},\hat{\bm{\lambda}}_{j})\end{bmatrix}+\begin{bmatrix}\mbox{\boldmath$\mathrm{K}$}_{uu}&\mbox{\boldmath$\mathrm{K}$}_{u\lambda}\\\
\mbox{\boldmath$\mathrm{K}$}_{\lambda
u}&\mbox{\boldmath$\mathrm{K}$}_{\lambda\lambda}\end{bmatrix}\begin{bmatrix}\Delta\hat{\bm{u}}\\\
\Delta\hat{\bm{\lambda}}\end{bmatrix}=\mbox{\boldmath$\mathrm{0}$}$

where the $\mathrm{K}$ terms abbreviate the current tangent stiffness matrix.
Their entries are calculated as

$\displaystyle\mbox{\boldmath$\mathrm{K}$}_{uu}=\frac{\partial\mbox{\boldmath$\mathrm{R}$}_{u}}{\partial\mbox{\boldmath$\mathrm{u}$}}=$
$\displaystyle\int_{\Omega_{e}}{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}\bigg{(}2\mu_{0}{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}({\mbox{\boldmath$\mathrm{F}$}^{e}}^{T}\mbox{\boldmath$\mathrm{F}$}^{e}-\mbox{\boldmath$\mathrm{C^{*}}$})$
$\displaystyle+4\mu_{0}{\mbox{\boldmath$\mathrm{F}$}^{e}}^{T}{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}\mbox{\boldmath$\mathrm{F}$}^{e}-\mu_{0}{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}\hat{\bm{\lambda}}{\mbox{\boldmath$\mathrm{B}$}^{e}}\hat{\bm{\lambda}}{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}\bigg{)}\,\mathrm{d}\Omega$
$\displaystyle\mbox{\boldmath$\mathrm{K}$}_{u\lambda}=\frac{\partial\mbox{\boldmath$\mathrm{R}$}_{u}}{\partial\bm{\lambda}}=$
$\displaystyle\int_{\Omega_{e}}-{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}\mbox{\boldmath$\mathrm{S}$}^{*}-{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}\mu_{0}{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}{\mbox{\boldmath$\mathrm{B}$}^{e}}\hat{\bm{\lambda}}{\mbox{\boldmath$\mathrm{F}$}^{e}}^{T}$
$\displaystyle-{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}\mu_{0}{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}\hat{\bm{\lambda}}{\mbox{\boldmath$\mathrm{B}$}^{e}}{\mbox{\boldmath$\mathrm{F}$}^{e}}^{T}\,\mathrm{d}\Omega$
$\displaystyle\mbox{\boldmath$\mathrm{K}$}_{\lambda
u}=\frac{\partial\mbox{\boldmath$\mathrm{R}$}_{\lambda}}{\partial\mbox{\boldmath$\mathrm{u}$}}=$
$\displaystyle\int_{\Omega_{e}}\mu_{0}{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}{\mbox{\boldmath$\mathrm{B}$}^{e}}{\mbox{\boldmath$\mathrm{B}$}^{e}}\hat{\bm{\lambda}}{\mbox{\boldmath$\mathrm{F}$}^{e}}^{T}+\mu_{0}{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}{\mbox{\boldmath$\mathrm{F}$}^{e}}{\mbox{\boldmath$\mathrm{B}$}^{e}}\hat{\bm{\lambda}}{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}$
$\displaystyle+{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}\mbox{\boldmath$\mathrm{S}$}^{*}\,\mathrm{d}\Omega$
$\displaystyle\mbox{\boldmath$\mathrm{K}$}_{\lambda\lambda}=\frac{\partial\mbox{\boldmath$\mathrm{R}$}_{\lambda}}{\partial\bm{\lambda}}=$
$\displaystyle\int_{\Omega_{e}}\mu_{0}{\mbox{\boldmath$\mathrm{B}$}^{e}}^{T}{\mbox{\boldmath$\mathrm{F}$}^{e}}{\mbox{\boldmath$\mathrm{B}$}^{e}}{\mbox{\boldmath$\mathrm{F}$}^{e}}^{T}\,\mathrm{d}\Omega$

In each time or load step and in each data iteration the iterative solution of
the system is necessary.

###  A.3 Data generation of the 3-D example

For the three dimensional example in Section 4.3 the six deformations of case
C are described by six homogenized strain values. Remind that in the linear
regime $\bm{C}\approx 2\bm{\epsilon}+1$ and we can also gain data of other
strain and stress equivalents in this regime. Therefore, we adopt small strain
notation which leads to the six strain tensors in Voigt notation

$\displaystyle\bar{\bm{\epsilon}}^{(1)}=\begin{pmatrix}\alpha\\\ 0\\\ 0\\\
0\\\ 0\\\ 0\end{pmatrix},\bar{\bm{\epsilon}}^{(2)}=\begin{pmatrix}0\\\
\alpha\\\ 0\\\ 0\\\ 0\\\
0\end{pmatrix},\bar{\bm{\epsilon}}^{(3)}=\begin{pmatrix}0\\\ 0\\\ \alpha\\\
0\\\ 0\\\ 0\end{pmatrix},\bar{\bm{\epsilon}}^{(4)}=\begin{pmatrix}0\\\ 0\\\
0\\\ \alpha\\\ 0\\\
0\end{pmatrix},\bar{\bm{\epsilon}}^{(5)}=\begin{pmatrix}0\\\ 0\\\ 0\\\ 0\\\
\alpha\\\ 0\end{pmatrix},\bar{\bm{\epsilon}}^{(6)}=\begin{pmatrix}0\\\ 0\\\
0\\\ 0\\\ 0\\\ \alpha\end{pmatrix}\,.$ (27)

with an arbitrary deformation defined by $\alpha\in\mathbb{R}$. The six
computations with the introduced RVE lead with $\alpha=0.02$ to the six
homogenized stresses $\sigma_{i}$ $i=1,\dots,6$:

$\displaystyle\bar{\sigma}^{(1)}$
$\displaystyle=\begin{pmatrix}9.13&3.87&3.96&-0.00&0.00&-0.00\end{pmatrix}^{T}\cdot
10^{4}\text{ Pa}$ $\displaystyle\bar{\sigma}^{(2)}$
$\displaystyle=\begin{pmatrix}3.87&8.57&3.89&0.00&0.00&-0.00\end{pmatrix}^{T}\cdot
10^{4}\text{ Pa}$ $\displaystyle\bar{\sigma}^{(3)}$
$\displaystyle=\begin{pmatrix}3.96&3.89&8.89&-0.00&0.00&-0.00\end{pmatrix}^{T}\cdot
10^{4}\text{ Pa}$ $\displaystyle\bar{\sigma}^{(4)}$
$\displaystyle=\begin{pmatrix}-0.00&0.00&0.00&1.94&-0.00&0.00\end{pmatrix}^{T}\cdot
10^{4}\text{ Pa}$ $\displaystyle\bar{\sigma}^{(5)}$
$\displaystyle=\begin{pmatrix}0.00&0.00&0.00&0.00&1.96&0.00\end{pmatrix}^{T}\cdot
10^{4}\text{ Pa}$ $\displaystyle\bar{\sigma}^{(6)}$
$\displaystyle=\begin{pmatrix}-0.00&-0.00&-0.00&0.00&-0.00&1.99\end{pmatrix}^{T}\cdot
10^{4}\text{ Pa}$