Finite Strain Homogenization Using a Reduced Basis and Efficient Sampling

Abstract

The computational homogenization of hyperelastic solids in the geometrically nonlinear context has yet to be treated with sufficient efficiency in order to allow for real-world applications in true multiscale settings. This problem is addressed by a problem-specific surrogate model founded on a reduced basis approximation of the deformation gradient on the microscale. The setup phase is based upon a snapshot POD on deformation gradient fluctuations, in contrast to the widespread displacement-based approach. In order to reduce the computational offline costs, the space of relevant macroscopic stretch tensors is sampled efficiently by employing the Hencky strain. Numerical results show speed-up factors in the order of 5–100 and significantly improved robustness while retaining good accuracy. An open-source demonstrator tool with 50 lines of code emphasizes the simplicity and efficiency of the method.

1. Introduction

1.1. Purpose

The description of solid mechanics under finite strains is of particular interest in both academia and industry. It allows for accurate descriptions of rotations and stretches under mild assumptions. Thus, many geometric effects can be captured. For instance, alignments and rearrangements of the respective structures may trigger pronounced stiffening or softening effects.

In such cases where rotations and deformations are not suitable for linearization, dissipative effects also play a notable role for many materials. Regardless of the kind of dissipation involved in a certain process, hyperelasticity usually persists to a certain extent. Therefore, it is worthwhile investigating this comparatively simple case at first, before introducing history dependence into the description. Prominent examples of materials that require a hyperelastic description at finite strains include carbon black-filled rubber [1] and amorphous glassy polymers [2], to name just two.

The main purpose of this work is the computationally efficient quasi-static homogenization of hyperelastic solids with full account for geometric nonlinearities. The employed methodology is twofold. First, a Reduced Basis (RB) model for the microscopic problem is established. The term Reduced Basis, used in this work, is not to be confused with the homonymous method introduced by Barrault, Maday, Nguyen, and Patera [3]. Once set up, it enables more efficient evaluations of the homogenized material response as compared to the Finite Element Method (FEM). Second, an efficient strategy for sampling of the space of macroscopic kinematic states is proposed. This renders the setup phase of the RB model more rational.

1.2. State of the Art

Efficiently determining the overall solid–mechanical properties of microstructures has been investigated for decades, and a large body of literature is available. Comprehensive review articles, such as [4] and [5], summarize the progress. Here, attention is restrained to few methods most similar or relevant to the present work.

The FE² method [6] is theoretically capable of performing realistic two-scale simulations with arbitrary accuracy. Therefore, it serves as a reference method in the context of first-order homogenization based on the assumption of separated length scales. In the FE², the evaluations of the unknown macroscopic constitutive law are approximated by microscopic FE simulations. However, this comes along with computational costs that quickly exceed the capabilities of common workstations, both at present and in the foreseeable future. Roughly speaking, the computational effort required on the microscale multiplies with that of the macroscale, hence the method’s name. It is thus worthwhile to develop order reduction methods for the microscopic problem.

A common approach within the field of computational homogenization (and well beyond) is to extract essential information from provided in silico data. To this end, schemes based on the Proper Orthogonal Decomposition (POD) compute correlations within snapshot data, [7]. Such methods include the R3M [8] and can be further enhanced by the use of, e.g., the EIM, as in [9]. Numerical comparisons of various schemes were conducted in [10,11]. To the best of the authors’ knowledge, all published POD-based methods addressing the finite strain hyperelastic problem choose to reduce the number of degrees of freedom (DOF) of the displacement field. This results in sometimes significant speed-ups. Another important feature is that they allow for reconstruction of the microscopic displacement fields. The application of the snapshot POD to gradients of the primal variables has been studied—e.g., for infinitesimal strain hyperelasticity [12] and fluid mechanics [13]—but does not appear to have been investigated for finite strain hyperelasticity yet.

Still, the solution of the reduced equations remains a complex task. It requires evaluations of material laws and numerical integration over the microstructure. Promising progress has been made in the field of efficient integration schemes, see for instance [14,15]. A main reason for the speed-up of these methods is the reduced number of function evaluations.

The highest speed-ups are achievable if the computational effort of the determination of effective microstructural responses can be fully decoupled from underlying microstructural discretizations. Such homogenization methods directly approximate the effective material law by means of a dedicated numerical scheme. Technically, this can be seen as the direct surrogation of unknown functions, e.g., of the effective free energy or stress. For instance, the Material Map [16] interpolates the coefficients of an assumed macroscopic material model. Another example is the NEXP method [17], where the effective stored energy density is approximated using a tensor product of one-dimensional splines. The authors treated the case of small strains by introducing the RNEXP method [12], where the effective stored energy is interpolated by a dedicated kernel scheme.

However, interpolatory and regressional methods suffer the inherent drawback of not providing any explicit information on the microscale. For instance, microscopic displacement or stress fields cannot be reconstructed from the solutions of macroscopic interpolation. Another important open question is how to provide the supporting data points for the interpolation in an efficient manner. The data at these points is usually provided by the solution of a full-order model (FOM) and come along with the corresponding numerical costs. Hence, the positions of data points in the parameter space should be chosen carefully, as unnecessary or redundant solutions of the FOM should be avoided. On the other hand, too sparsely seeded points might not capture the homogenized properties of the microstructure appropriately.

1.3. Main Contributions and Outline

The present work generalizes parts of the previous paper [12] to the finite strain regime. It aims at reducing the computational complexity for the determination of the homogenized microstructural response, which is parametrized by the macroscopic deformation gradient acting as a boundary condition. This is achieved by means of a Reduced Basis approximation of the microscopic deformation gradient. The basis is obtained with the aid of a POD of snapshots of fluctuation fields of the deformation gradient. Thus, the application of the RB model does not necessitate the computation of gradients of displacement fields, and even does not require the displacements to be available at all. In other words, microscopic displacement fields are completely avoided. However, they can be reconstructed from the RB approximation of the deformation gradient, uniquely up to rigid body motion.

Another key advantage is the sleek implementation of the method. A demonstration containing a minimum working example of the RB model with 50 lines of MATLAB/Octave code is provided [18].

As for the setup phase, the snapshot data is created by means of an efficient sampling procedure for the microscopic boundary condition. To this end, the set of macroscopic Hencky strains is identified as a suitable linear parameter space, within which the sampling sites are placed based upon physical interpretation. This allows for control of the resolution of certain key characteristics of the effective material response while keeping the total number of samples within bounds.

The Reduced Basis method is presented in Section 2. The basis identification is based on the sampling strategy developed in Section 3. Numerical examples are presented in Section 4. Both the numerical and the theoretical findings are summarized and discussed in Section 5.

1.4. Notation

The set of real numbers and the subset of positive numbers greater than zero are denoted by $\mathbb{R}$ and ${\mathbb{R}}_{+}$, respectively. Matrices are marked by two underlines and vectors by one underline, e.g., $\underline{\underline{A}}$, $\underline{a}$. Vectors are assumed to be columns, and the dot product of two vectors of the same size is understood as the Euclidean scalar product, $\underline{x}·\underline{y}={\underline{x}}^{\mathrm{T}}\underline{y}$. First order and second order tensors in coordinate-free description are denoted by bold letters, e.g., A, a. No conclusion can be drawn on the order of a tensor based on its capitalization. Here, the underlying space is always the Euclidean space ${\mathbb{R}}^3$ with its standard basis. First order and second order tensors can be represented as vectors and matrices, e.g., $\mathit{A}\leftrightarrow \underline{A}\in {\mathbb{R}}^3$ and $\mathit{B}\leftrightarrow \underline{\underline{B}}\in {\mathbb{R}}^{3\times 3}$, respectively. Norms of vectors and matrices respectively denote the Euclidean and the Frobenius norm. The norm of a tensor of second order equals the norm of its matrix representation for the chosen basis. Fourth order tensors are denoted by blackboard bold symbols other than $\mathbb{R}$, e.g., $\mathbb{C}$ and $\mathbb{I}$. Components of tensors of order M are with respect to the Euclidean tensorial basis ${\mathit{e}}^{(1)}\otimes {\mathit{e}}^{(2)}\otimes \cdots \otimes {\mathit{e}}^{(M)}$, e.g., $A_{ij}$, $B_{ij}$ for second order tensors A, B and $C_{ijkl}$, $C_{ijkl}^{\prime }$ for $\mathbb{C}$, ${\mathbb{C}}^{\prime }$. The following contractions are defined:

$$\begin{array}{cccc}\hfill \mathit{A}·\mathit{B}& =\sum _{i,j=1}^3{A}_{ij}{B}_{ij},\hfill & \hfill \mathbb{C}·\mathit{B}& =\sum _{i,j,k,l=1}^3{C}_{ijkl}{B}_{kl}{\mathit{e}}^{(i)}\otimes {\mathit{e}}^{(j)},\hfill \\ \hfill \mathit{A}·\mathbb{C}& =\sum _{i,j,k,l=1}^3{A}_{ij}{C}_{ijkl}{\mathit{e}}^{(k)}\otimes {\mathit{e}}^{(l)},\hfill & \hfill \mathbb{C}·{\mathbb{C}}^{\prime }& =\sum _{i,j,k,l,m,n=1}^3{C}_{ijmn}{C}_{mnkl}^{\prime }{\mathit{e}}^{(i)}\otimes {\mathit{e}}^{(j)}\otimes {\mathit{e}}^{(k)}\otimes {\mathit{e}}^{(l)}.\hfill \end{array}$$

Let $\Omega \subset {\mathbb{R}}^3$ be the domain occupied by a physical body undergoing elastic deformations, and let ${\Omega }_0$ be its initial configuration. Then, x and X describe the coordinates of material points within the current configuration $\Omega$ and within the reference state ${\Omega }_0$, respectively. Their difference is the displacement $\mathit{u}=\mathit{x}-\mathit{X}$, see Figure 1. The gradient of a vector field $v=v(\mathit{X})$ is defined as a right gradient and denoted by ${\displaystyle \frac{{\displaystyle \partial v}}{{\displaystyle \partial \mathit{X}}}=v\otimes {\mathsf{\nabla }}_{\mathrm{X}}}$. The divergence of a second order tensor field is the vector field resulting from row-wise divergence. The boundaries of the respective configurations are denoted by $\partial \Omega$ and $\partial {\Omega }_0$. The set of square-integrable Lebesgue functions on the reference domain is tagged $L^2({\Omega }_0)$.

The displacement gradient $\mathit{H}=\mathit{u}\otimes {\mathsf{\nabla }}_{\mathrm{X}}$ and the deformation gradient $\mathit{F}=\mathit{x}\otimes {\mathsf{\nabla }}_{\mathrm{X}}$ are related through $\mathit{F}=\mathit{H}+\mathit{I}$, where I is the second order identity tensor in three dimensions. The determinant $J=det\mathit{F}$ measures the relative volumetric change due to the present deformation.

Unimodular quantities, i.e., second order tensors with determinant ones, may be emphasized by a hat, e.g., $\widehat{\mathit{F}}=J^{-1/3}\mathit{F}$. This multiplicative decomposition is sometimes attributed to Flory [19] and also goes by the name Dilatational-Deviatoric Multiplicative Split (DDMS).

In the two-scale context, overlined symbols represent quantities on the macroscopic scale, e.g., $\overline{\mathit{A}}$, $\overline{\mathit{a}}$, while symbols without overline correspond to their microscopic counterpart, e.g., A, a. Equivalently, macroscopic quantities are called global and microscopic ones are called local. The volume average of a general local field $\phi$

$$\begin{array}{cc}\hfill 〈\phi 〉& =〈\phi (•)〉=\frac{1}{|{\Omega }_0|}{\int }_{{\Omega }_0}\phi (•)\mathrm{d}V\hfill \end{array}$$

(1)

is essential to the theory. The dependence of a microscopic quantity A on both the microscopic coordinates X and a macroscopic quantity $\overline{\mathit{B}}$ is denoted by $\mathit{A}=\mathit{A}(\mathit{X};\overline{\mathit{B}})$. In such a case, the components of the macroscopic quantity $\overline{\mathit{B}}$ are called parameters of the microscopic function $\mathit{A}(•;\overline{\mathit{B}})$. The application of the volume averaging operator is abbreviated by $〈\mathit{A}〉=〈\mathit{A}(•;\overline{\mathit{B}})〉$. The case of a concatenated function $f(\mathit{A})=f(\mathit{A}(\mathit{X};\overline{\mathit{B}}))$ is analogous, i.e. $〈f〉=〈f(\mathit{A})〉=〈f(\mathit{A}(•;\overline{\mathit{B}}))〉$, regardless of the tensorial order of the image of the function f.

1.5. Material Models

In this work, hyperelastic materials are investigated. They are characterized by stored energy density functions $W=W(\mathit{F})$. The first Piola–Kirchhoff stress

$$\begin{array}{cc}\hfill \mathit{P}(\mathit{F})& {\displaystyle =\frac{{\displaystyle \partial W}}{{\displaystyle \partial \mathit{F}}}(\mathit{F})}\hfill \end{array}$$

(2)

and the corresponding fourth-order stiffness tensor

$$\begin{array}{cc}\hfill \mathbb{C}(\mathit{F})& {\displaystyle =\frac{{\displaystyle {\partial }^{2}W}}{{\displaystyle \partial {\mathit{F}}^2}}(\mathit{F})}\hfill \end{array}$$

(3)

characterize the material response.

Henceforth, for reasons of readability, the stored energy density function W will be spoken of as an energy, and the terms stored and density will not always be mentioned explicitly. In the infinitesimal strain framework, hyperelastic energies have been formulated to model deformation plasticity (e.g., [12,17,20]). Although these models are only valid for purely proportional loading conditions, they provide means to simulate highly nonlinear material behavior in certain scenarios comparably easily within the context of hyperelasticity. Note that genuine dissipative processes require additional state describing variables with corresponding evolution laws.

The proposed method is suitable for any type of hyperelastic constitutive law. As the modeling of complex material behavior is not the main focus of this study, the Neo–Hookean law

$$\begin{array}{cc}\hfill W(\mathit{F})=W_{\mathrm{DDMS}}(J,\widehat{\mathit{F}})& =\frac{K}{4}\left[{(J-1)}^2+{(\ln J)}^2\right]+\frac{G}{2}\left(\mathrm{tr}({\widehat{\mathit{F}}}^{\mathsf{T}}\widehat{\mathit{F}})-3\right)\hfill \end{array}$$

(4)

is used, with K the bulk modulus and G the shear modulus. The volumetric part of the energy is taken from [21]. Using the DDMS, a decoupled dependence on the volumetric and isochoric part of the deformation is assumed, which is a common way to model the distinct material behavior with respect to these two contributions, see e.g., [22].

1.6. Problem Setting of First Order Homogenization

Assuming stationarity and separability of scales, the following coupled and deformation-driven problems can be derived by means of asymptotic expansion of the displacement u and subsequent first order approximation. This procedure is carried out in [23] with much detail. Here, the technical process is omitted and only the resulting equations are stated.

1.6.1. Macroscopic Problem

Balance of linear momentum

$$\begin{array}{cc}\hfill {\mathrm{Div}}_{\overline{\mathrm{X}}}(\overline{\mathit{P}})+\overline{\mathit{b}}& =0,\hfill \end{array}$$

(5)

where $\overline{\mathit{b}}$ denote bulk forces, and balance of angular momentum

$$\begin{array}{cc}\hfill {\overline{\mathit{F}}}^{-1}\overline{\mathit{P}}& ={\overline{\mathit{P}}}^{\mathsf{T}}{\overline{\mathit{F}}}^{-\mathsf{T}},\hfill \end{array}$$

(6)

along with well-posed boundary conditions that constitute the macroscopic boundary value problem. This system of equations is closed by means of the hyperelasticity law, cf. (2),

$$\begin{array}{cc}\hfill \overline{\mathit{P}}(\overline{\mathit{F}})& {\displaystyle =\frac{{\displaystyle \partial \overline{W}}}{{\displaystyle \partial \overline{\mathit{F}}}}.}\hfill \end{array}$$

(7)

Note that, in general, $\overline{W}$ is a priori not available and $(7)$ is thus a purely formal relation. For reasons of readability, the dependence of any quantity on the macroscopic material coordinate $\overline{\mathit{X}}$ is usually spared, e.g., $\overline{\mathit{F}}=\overline{\mathit{F}}(\overline{\mathit{X}})$, $\overline{\mathit{H}}=\overline{\mathit{H}}(\overline{\mathit{X}})$, or $\overline{\mathit{u}}=\overline{\mathit{u}}(\overline{\mathit{X}})$.

1.6.2. Microscopic Problem

The microscopic boundary value problem is given by the balance equations

$$\begin{array}{cc}\hfill {\mathrm{Div}}_{\mathrm{X}}(\mathit{P})& =0\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill {\mathit{F}}^{-1}\mathit{P}& ={\mathit{P}}^{\mathsf{T}}{\mathit{F}}^{-\mathsf{T}}\hfill \end{array}$$

(9)

and suitable boundary conditions, e.g., as discussed in [24]. In this work, periodic displacement fluctuation boundary conditions are employed. The microscopic displacements take the form

$$\begin{array}{cc}\hfill \mathit{u}(\mathit{X};\overline{\mathit{F}})& =\overline{\mathit{u}}+\overline{\mathit{H}}\mathit{X}+\mathit{w}(\mathit{X};\overline{\mathit{F}}).\hfill \end{array}$$

(10)

Therein, the macroscopic displacement is independent of microscopic quantities. The second term, $\overline{\mathit{H}}\mathit{X}$, corresponds to a homogeneous deformation of the microstructure. The third term, $\mathit{w}(\mathit{X})$, is a displacement fluctuation with the zero mean property $\langle \mathit{w}\rangle =0$. Hence, the deformation gradient reads

$$\begin{array}{cc}\hfill \mathit{F}(\mathit{X};\overline{\mathit{F}})& =\overline{\mathit{F}}+\tilde{\mathit{H}}(\mathit{X};\overline{\mathit{F}})=\overline{\mathit{F}}+\tilde{\mathit{F}}(\mathit{X};\overline{\mathit{F}}),\hfill \end{array}$$

(11)

where the fluctuation part $\tilde{\mathit{H}}=\mathit{w}\otimes {\mathsf{\nabla }}_{\mathrm{X}}=\tilde{\mathit{F}}$, has the zero mean property

$$\begin{array}{cc}\hfill 〈\tilde{\mathit{F}}〉& =0.\hfill \end{array}$$

(12)

Thus, the volume average of the local deformation gradient equals its macroscopic counterpart,

$$\begin{array}{cc}\hfill \overline{\mathit{F}}& =〈\mathit{F}〉.\hfill \end{array}$$

(13)

This motivates the identification of $\overline{\mathit{F}}$ as the boundary condition to the microscopic problem (8). As for the material response, the following relationships can be deduced:

$$\begin{array}{cc}\hfill \overline{W}& =〈W〉,\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill \overline{\mathit{P}}& =〈\mathit{P}〉,\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill \overline{\mathbb{C}}& \ne 〈\mathbb{C}〉.\hfill \end{array}$$

(16)

Equations (13) and (15) are called kinematic and static coupling relations, respectively. The inequality (16) generally holds for heterogeneous microstructures, even in the most simple case of infinitesimal strains and linear elasticity. More precisely, the volume average overestimates the effective stiffness in the spectral sense,

$$\begin{array}{cc}\hfill \mathit{x}·〈\mathbb{C}〉·\mathit{x}& \ge \mathit{x}·\overline{\mathbb{C}}·\mathit{x}\forall \mathrm{second}\mathrm{order}\mathrm{tensors}\mathit{x}.\hfill \end{array}$$

(17)

2. Reduced Basis Homogenization for Hyperelasticity

2.1. Formulation

The Reduced Basis (RB) scheme is based on a direct approximation of the microscopic deformation gradient F from Equation (11) without the need to explicitly have the corresponding displacement at hand. The initial approximation is given by

$$\begin{array}{cc}\hfill {\mathit{F}}_{\xi }(\mathit{X};\overline{\mathit{F}},\underline{\xi })& =\overline{\mathit{F}}+\sum _{i=1}^N{\mathit{B}}^{(i)}(\mathit{X}){\xi }_i.\hfill \end{array}$$

(18)

The arguments $\overline{\mathit{F}}$ and $\underline{\xi }\in {\mathbb{R}}^N$ are the boundary condition to (8) and the reduced coefficient vector, respectively. Note that the macroscopic coordinate $\overline{\mathit{X}}$ is not assumed to influence the RB approximation, i.e., the same approximation is made throughout the macrostructure. The set ${\{{\mathit{B}}^{(i)}\}}_{i=1}^N$ is linearly independent within the space $L^2({\Omega }_0)$ and is called RB of F. In a later context, it will also be referred to as the set of ansatz functions. In order to enforce the relationship

$$\begin{array}{cc}\hfill 〈{\mathit{F}}_{\xi }〉& =\overline{\mathit{F}}\hfill \end{array}$$

(19)

regardless of the reduced coefficients $\underline{\xi }$, the basis functions are asserted to have the fluctuation property, i.e., for $i=1,\dots ,N$

$$\begin{array}{cc}\hfill 〈{\mathit{B}}^{(i)}〉& =0.\hfill \end{array}$$

(20)

For now, the RB is assumed to be given.

The ansatz (18) allows for evaluation of the energy at the macroscale as a function of the macroscopic kinematic variable $\overline{\mathit{F}}$ and of the reduced coefficients $\underline{\xi }$,

$$\begin{array}{cc}\hfill {\overline{W}}_{\xi }(\overline{\mathit{F}},\underline{\xi })& =〈W({\mathit{F}}_{\xi })〉.\hfill \end{array}$$

(21)

By the principle of minimum energy, the optimal coefficients

$$\begin{array}{cc}\hfill {\underline{\xi }}^{*}(\overline{\mathit{F}})& =arg\underset{\underline{\xi }\in {\mathbb{R}}^N}{min}{\overline{W}}_{\xi }(\overline{\mathit{F}},\underline{\xi })\hfill \end{array}$$

(22)

are sought after. The unconstrained and unique solvability of this task is assumed for the moment and will be addressed in Section 2.4. The solution of (22) defines the RB approximation of the deformation gradient

$$\begin{array}{cc}\hfill {\mathit{F}}^{\mathrm{RB}}(\mathit{X};\overline{\mathit{F}})& =\overline{\mathit{F}}+\sum _{i=1}^N{\mathit{B}}^{(i)}(\mathit{X}){\xi }_i^{*}(\overline{\mathit{F}}).\hfill \end{array}$$

(23)

The microscopic energy, stress, and stiffness within the microstructure are then approximated by

$$\begin{array}{cc}\hfill W^{\mathrm{RB}}(\mathit{X};\overline{\mathit{F}})& =W({\mathit{F}}^{\mathrm{RB}}(\mathit{X};\overline{\mathit{F}})),\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill {\mathit{P}}^{\mathrm{RB}}(\mathit{X};\overline{\mathit{F}})& {\displaystyle =\frac{{\displaystyle \partial W}}{{\displaystyle \partial \mathit{F}}}({\mathit{F}}^{\mathrm{RB}}(\mathit{X};\overline{\mathit{F}})),}\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill {\mathbb{C}}^{\mathrm{RB}}(\mathit{X};\overline{\mathit{F}})& {\displaystyle =\frac{{\displaystyle {\partial }^{2}W}}{{\displaystyle \partial {\mathit{F}}^2}}({\mathit{F}}^{\mathrm{RB}}(\mathit{X};\overline{\mathit{F}})),}\hfill \end{array}$$

(26)

respectively. The resulting effective energy is readily given by

$$\begin{array}{cc}\hfill {\overline{W}}^{\mathrm{RB}}(\overline{\mathit{F}})& =〈W^{\mathrm{RB}}〉.\hfill \end{array}$$

(27)

Also, the effective responses ${\overline{\mathit{P}}}^{\mathrm{RB}}(\overline{\mathit{F}})$ and ${\overline{\mathbb{C}}}^{\mathrm{RB}}(\overline{\mathit{F}})$ may now be calculated. However, before going into detail on that, it is advantageous to first elaborate on the solution of the minimization problem (22). This short survey will reveal essential properties of some occurring quantities that are important for the determination of the effective material response.

The necessary, first order optimality conditions to (22) define the components of the residual vector $\underline{r}\in {\mathbb{R}}^N$,

$$\begin{array}{cc}\hfill r_i(\overline{\mathit{F}},\underline{\xi })& {\displaystyle =\frac{{\displaystyle \partial {\overline{W}}_{\xi }}}{{\displaystyle \partial {\xi }_i}}(\overline{\mathit{F}},\underline{\xi })=〈{\displaystyle \frac{{\displaystyle \partial W}}{{\displaystyle \partial \mathit{F}}}({\mathit{F}}_{\xi })·\frac{{\displaystyle \partial {\mathit{F}}_{\xi }}}{{\displaystyle \partial {\xi }_i}}}〉=〈\mathit{P}({\mathit{F}}_{\xi })·{\mathit{B}}^{(i)}〉=0(i=1,\dots ,N).}\hfill \end{array}$$

(28)

Note 1.

The solution stress field ${\mathit{P}}^{\mathrm{RB}}$ is $L^2({\Omega }_0)$-orthogonal to the RB ansatz functions ${\{{\mathit{B}}^{(i)}\}}_{i=1}^N$.

A viable choice for the solution of the minimization problem (22) is the Newton–Raphson scheme, which necessitates the Jacobian $\underline{\underline{D}}\in \mathrm{Sym}({\mathbb{R}}^{N\times N})$ with the components

$$\begin{array}{cc}\hfill D_{ij}& =〈{\displaystyle {\mathit{B}}^{(i)}·\frac{{\displaystyle \partial \mathit{P}}}{{\displaystyle \partial {\xi }_j}}({\mathit{F}}_{\xi })}〉=〈{\mathit{B}}^{(i)}·\mathbb{C}({\mathit{F}}_{\xi })·{\mathit{B}}^{(j)}〉=D_{ji}(i,j=1,\dots ,N).\hfill \end{array}$$

(29)

Then, the kth iteration to the solution ${\underline{\xi }}^{*}(\overline{\mathit{F}})$ reads

$$\begin{array}{cc}\hfill {\underline{\xi }}^{(k)}& ={\underline{\xi }}^{(k-1)}-{\left({\underline{\underline{D}}}^{(k-1)}\right)}^{-1}{\underline{r}}^{(k-1)}(k\ge 1).\hfill \end{array}$$

(30)

The initial guess ${\underline{\xi }}^{(0)}$ can be zero or a more sophisticated alternative.

The deduction of the effective material response by means of Note 1 and the Jacobian $\underline{\underline{D}}$ is given in Appendix B. The following list summarizes the homogenized quantities arising from the F-RB:

$$\begin{array}{ccc}\hfill \overline{\mathit{F}}& =〈{\mathit{F}}^{\mathrm{RB}}〉\hfill & \mathrm{see}(19)\hfill \\ \hfill {\overline{W}}^{\mathrm{RB}}& =〈W^{\mathrm{RB}}〉\hfill & \mathrm{see}(27)\hfill \\ \hfill {\overline{\mathit{P}}}^{\mathrm{RB}}& =〈{\mathit{P}}^{\mathrm{RB}}〉\hfill & \mathrm{see}\mathrm{Appendix}\mathrm{B}.1\hfill \\ \hfill {\overline{\mathbb{C}}}^{\mathrm{RB}}& =〈{\mathbb{C}}^{\mathrm{RB}}〉-\sum _{i,j=1}^N{D}_{ij}^{-1}〈{\mathbb{C}}^{\mathrm{RB}}·{\mathit{B}}^{(i)}〉\otimes 〈{\mathit{B}}^{(j)}·{\mathbb{C}}^{\mathrm{RB}}〉\hfill & \mathrm{see}\mathrm{Appendix}\mathrm{B}.2\hfill \end{array}$$

In total, the problem of determining both the local field F and the homogenized material properties depends only on N degrees of freedom, namely on the N coefficients ${\xi }_i$. Usually, the number of DOF N is in the range of 10 to 150, which compares to the full order model’s complexity that can easily reach more than 10⁵ DOF.

Remark 1.

Despite this impressive reduction of the number of DOF, the computational costs associated with the assembly of the residual $\underline{r}$ and of the Jacobian $\underline{\underline{D}}$ still relate to the number of quadrature points of the microstructural discretization.

This method extends corresponding methods for the geometrically linear case, where the infinitesimal strain tensor $\mathit{\epsilon }=\mathrm{sym}(\mathit{H})$ is considered. For more information on these topics, the reader is referred to the authors’ previous work [12] or standard literature, such as ([25], part II.C).

2.2. Identification of the Reduced Basis

The basis ${\{{\mathit{B}}^{(i)}\}}_{i=1}^N$ is computed by means of a classical snapshot POD. In contrast to many other POD based reduction methods, it is important to point out that here, the primal variable is not taken to be the displacement field, u. Instead, the POD is performed on deformation gradient fluctuations, $\tilde{\mathit{F}}$.

During the snapshot POD, snapshots are first created by means of high-fidelity solutions to the nonlinear microscopic problem (8) with different snapshot boundary conditions ${\overline{\mathit{F}}}^{(j)}$, $j=1,\dots ,N_{\mathrm{s}}$, which are also called training boundary conditions. Each of these boundary conditions leads to a solution field ${\mathit{F}}^{(j)}(\mathit{X};\overline{\mathit{F}})$. Typically, the Finite Element Method (FEM) or solvers making use of the Fast Fourier Transform (e.g., [26]) are used to this end. The RB method presented here is independent of the discretization method utilized to obtain full field solutions. It is merely necessary to know the quadrature weights and the related discrete values of the solutions ${\mathit{F}}^{(j)}(\mathit{X};{\overline{\mathit{F}}}^{(j)})$. For better readability, the continuous notation is maintained for the moment. The corresponding fluctuation fields are computed by means of local subtraction of the macroscopic deformation gradient

$$\begin{array}{cc}\hfill {\tilde{\mathit{F}}}^{(j)}(\mathit{X};{\overline{\mathit{F}}}^{(j)})& ={\mathit{F}}^{(j)}(\mathit{X};{\overline{\mathit{F}}}^{(j)})-{\overline{\mathit{F}}}^{(j)}(j=1,\dots ,N_{\mathrm{s}}).\hfill \end{array}$$

(31)

Each of these $N_{\mathrm{s}}$ fluctuation fields ${\tilde{\mathit{F}}}^{(j)}$ represent one snapshot.

Second, the most dominant features of the snapshots are extracted. This is done by means of the eigendecomposition of the correlation matrix. It consists of the mutual $L^2({\Omega }_0)$ scalar products of the snapshots, $〈{\tilde{\mathit{F}}}^{(i)}·{\tilde{\mathit{F}}}^{(j)}〉$ ($i,j=1,\dots ,N_{\mathrm{s}}$), cf. (1). The remaining procedure is standard, see for instance [7] or [27]: the $N_{\mathrm{s}}$ eigenvalues ${\lambda }_j$, corresponding to the eigenvectors ${\underline{E}}^{(j)}\in {\mathbb{R}}^{N_{\mathrm{s}}}$, are sorted in descending order and truncated by only considering the first N values, ${\lambda }_1\ge \dots \ge {\lambda }_N$. The decision on a particular threshold index N is based on consideration of the Schmidt–Eckhard–Young–Mirsky theorem. Finally, the RB is constructed as

$$\begin{array}{cc}\hfill {\mathit{B}}^{(i)}(\mathit{X})& =\sum _{j=1}^{N_{\mathrm{s}}}\frac{1}{\sqrt{{\lambda }_i}}{\left({\underline{E}}^{(i)}\right)}_j{\tilde{\mathit{F}}}^{(j)}(\mathit{X})(i=1,\dots ,N)\hfill \end{array}$$

(32)

where the factor $1/\sqrt{{\lambda }_i}$ accounts for $L^2({\Omega }_0)$ normalization of the base elements.We conclude that the RB is a collection of $L^2({\Omega }_0)$ orthonormal $\tilde{\mathit{H}}$-like quantities.

2.3. Mathematical Motivation of the Reduced Basis Model

Next, the obtained deformation gradient field ${\mathit{F}}^{\mathrm{RB}}$ and the associated stress field ${\mathit{P}}^{\mathrm{RB}}$ are shown to weakly solve the original problem (8), ${\mathrm{Div}}_{\mathrm{X}}(\mathit{P})=0$.

Let $\delta \mathit{w}\in {\mathit{H}}_0^1({\Omega }_0)$ be an admissible test function, i.e., a once weakly differentiable periodic displacement fluctuation field, and let $\delta \mathit{H}=\delta \mathit{w}\otimes {\mathsf{\nabla }}_{\mathrm{X}}$ denote its gradient. Then, the well-known weak form of (8) is equivalent to the principle of virtual work,

$$\begin{array}{cc}\hfill \delta \overline{W}=& {\int }_{{\Omega }_0}\mathit{P}·\delta \mathit{H}\mathrm{d}V=0.\hfill \end{array}$$

(33)

The residual $\underline{r}$ from (28) coincides with the integral of the weak form, if the test function $\delta \mathit{w}$ is chosen suitably. As the basis elements ${\mathit{B}}^{(i)}$ are linear combinations of deformation gradient fluctuations ${\mathit{w}}^{(j)}\otimes {\mathsf{\nabla }}_{\mathrm{X}}$, cf. (32), the equivalence of (28) and (33) is obvious.

It follows that the Reduced Basis scheme is a Galerkin procedure, taking the displacement fields corresponding to the RB elements ${\mathit{B}}^{(i)}$ as both ansatz and test functions. Hence, the RB model is equivalent to the FEM, but with basis functions that are globally supported in ${\Omega }_0\backslash \partial {\Omega }_0$. In other words, the basis functions of the RB method span a subspace of dimension N within the high-dimensional space of FE basis functions. Although this property is shared with RB schemes based on POD of displacement snapshots, a notable difference is that this novel approach directly operates on fields entering the constitutive equations.

2.4. Details on the Coefficient Optimization

The coefficient optimization task (22) leads to a weak solution of the microscopic boundary value problem, as was just shown. Hence, the well-established theories on which the FEM is built, e.g., the calculus of variations, are applicable to the presented method just as much. This implies that the well-known issues with suitable convexity conditions and with existence and uniqueness of minimizers apply to the RB method, too. We focus on ad hoc numerical treatments of these issues. For more details on the theoretical part, the reader is referred to standard literature, such as [28], and recent developments in this matter, e.g., [29].

A constraint to the optimization problem is the physical condition

$$\begin{array}{cc}\hfill J& =\det ({\mathit{F}}_{\xi }(\mathit{X}))>0\forall \mathit{X}\in {\Omega }_0,\hfill \end{array}$$

(34)

which means that no singular ($J=0$) or self-penetrating ($J<0$) deformations shall occur. This reduces the set of admissible coefficients $\underline{\xi }$ to a subset of ${\mathbb{R}}^N$ that is nontrivial to access. The positiveness of the microscopic determinant of the deformation gradient introduces a constraint to the, thus far, unconstrained minimum problem (22), representing the weak form of (8) in the RB setting.

In case of a violation of the inequality (34), the implementation of the RB method is prone to failure as soon as the constitutive law (4) is evaluated. This occurs only in the context of volume averaging, i.e., for the assembly of the residual, the Jacobian, or the effective energy, stress, or stiffness. The numerical quadrature approximating the volume averaging operation is

$$\begin{array}{cc}\hfill 〈•〉& =\frac{1}{|{\Omega }_0|}{\int }_{{\Omega }_0}•(\mathit{X})\mathrm{d}V\approx \frac{1}{|{\Omega }_0|}\sum _{p=1}^{N_{\mathrm{qp}}}•({\mathit{X}}_p){\mathit{w}}_p.\hfill \end{array}$$

(35)

Here, $N_{\mathrm{qp}}$, ${\mathit{X}}_p$, and ${\mathit{w}}_p$ respectively denote the number of quadrature points, their positions, and the corresponding positive weights. Moreover, even if the inequality (34) is satisfied everywhere, the local field ${\mathit{F}}_{\xi }$ might possibly have some positive but overly small values of the determinant, $0<\det ({\mathit{F}}_{\xi }(\mathit{X}))\ll 1$, that are unphysical. In such a case, the energy optimization, cf. (22), would be dominated by these nearly singular points. Instead of allowing the optimization to focus almost exclusively on these exceptional points, we interpret unphysically small values of the determinant as limitations to the reliability of the RB method. On the other hand, large values $\det ({\mathit{F}}_{\xi }(\mathit{X}))\gg 1$ are not too detrimental to the functionality of the scheme, although such values are just as questionable.

Thus, the following weighted numerical quadrature rule is introduced:

$$\begin{array}{cc}\hfill 〈•〉& \approx \left(\sum _{p=1}^{N_{\mathrm{qp}}}•({\mathit{X}}_p)\varphi (J_p){\mathit{w}}_p\right)/\left(\sum _{p=1}^{N_{\mathrm{qp}}}\varphi (J_p){\mathit{w}}_p\right).\hfill \end{array}$$

(36)

Therein, the almost smooth cutoff function $\varphi :\mathbb{R}\to {\mathbb{R}}_{\ge 0}$ is empirically defined by

$$\begin{array}{cc}\hfill \varphi (J)& =\left\{\begin{array}{cc}1\hfill & \mathrm{if}J>0.6\hfill \\ 0.5\mathrm{erf}(30J-15)+0.5\hfill & \mathrm{if}0.4<J\le 0.6\hfill \\ 0\hfill & \mathrm{if}J\le 0.4\hfill \end{array}\right..\hfill \end{array}$$

(37)

which makes use of the well-known error function. The cutoff function is depicted in Figure 2. This reliability indicator could, in principle, be modified, e.g., the steepness, the smoothness, and its center could be adjusted. Thus, it should be regarded as an example only.

This weighted numerical quadrature rule is used henceforward for all numerical volume averaging operations. Its application will not be noted explicitly. However, the theoretical derivation of the RB method, as described in Section 2.1, is not affected by this, i.e., volume averaging operations remain exact as far as the theory is concerned. The two most important consequences of this numerical tweak are:

• The RB method is robust with respect to outlier values of the determinant. The modified quadrature rule extends the set of coefficient vectors $\underline{\xi }$ for which effective quantities can be computed, albeit approximately, to the whole space ${\mathbb{R}}^N$.
• The significance of local fields varies with the value of the cutoff function. When $\varphi$ attains values less than one, information is considered accordingly less reliable. In this sense, microscopic information is filtered based on a trust region for J defined by $\varphi$ can be seen as a reliability indicator.

3. Sampling

3.1. General Considerations

The proposed sampling strategy builds on the previous work [12], in which the authors proposed an analogous sampling procedure for the small strain setting. However, substantial modifications are required in order to account for the finite rotational part, $\overline{\mathit{R}}$, of the macroscopic deformation gradient, $\overline{\mathit{F}}$, and the nonlinearity of the volumetric part of the deformation, J, with respect to the local displacements, u. For the setup of the Reduced Basis model as described in Section 2.2, the space of macroscopic deformation gradients,

$$\begin{array}{cc}\hfill \overline{\mathcal{F}}& =\{\mathrm{second}\mathrm{order}\mathrm{tensors}\overline{\mathit{F}}|det\overline{\mathit{F}}>0\},\hfill \end{array}$$

(38)

needs to be sampled, i.e., the discrete sampling set ${\overline{\mathcal{F}}}_{\mathrm{s}}=\{\overline{\mathit{F}}{\}}_{m=1}^{N_{\mathrm{s}}}\subset \overline{\mathcal{F}}$ is to be defined. Two contradictory requirements need to be satisfied when constructing ${\overline{\mathcal{F}}}_{\mathrm{s}}$:

• The samples should be densely and homogeneously distributed within the space of all admissible macroscopic kinematic configurations. This is owing to the desire that the POD may extract correlation information from a holistic and unbiased set. In other words, the samples should be as uniformly random as possible within the anticipated query domain of the surrogate.
• The sample number $N_{\mathrm{s}}$ should not exceed a certain limit. Only with this property may the RB be identified within the bounds of available computational resources (e.g., memory and CPU time).

3.2. Large Strain Sampling Strategy

The set of admissible macroscopic deformation gradients $\overline{\mathcal{F}}$ is a subset of a nine-dimensional space ($\overline{\underline{\underline{\mathit{F}}}}\in {\mathbb{R}}^{3\times 3}\sim {\mathbb{R}}^9$) restricted by the inequality

$$\begin{array}{c}\hfill det\overline{\mathit{F}}>0.\end{array}$$

(39)

Therefore, a regular grid in the components of $\overline{\mathit{F}}$ might lead to a prohibitively large amount of samples, and even to a violation of (39). For instance, such a grid with a rather moderate resolution of just 10 samples of each component would require 1 billion solutions of the FOM. Also, the subsequent POD would involve a snapshot matrix and/or correlation matrix of accordingly huge dimensionality.

In order to decrease the dimension of the sampling space, recall the polar decomposition of the deformation gradient, $\overline{\mathit{F}}=\overline{\mathit{R}}\overline{U}$, where $\overline{\mathit{R}}$ is a rotation and $\overline{\mathit{U}}$ is the symmetric positive definite (s.p.d.) stretch tensor. Material objectivity implies the energy function to be independent of the frame of reference,

$$\begin{array}{cc}\hfill \overline{W}(\overline{\mathit{R}}\overline{\mathit{U}})& =\overline{W}(\overline{\mathit{U}}),\hfill \end{array}$$

(40)

and the transformation behavior

$$\begin{array}{cc}\hfill \overline{\mathit{P}}(\overline{\mathit{R}}\overline{\mathit{U}})& =\overline{\mathit{R}}\overline{\mathit{P}}(\overline{\mathit{U}})\hfill \end{array}$$

(41)

for the first Piola–Kirchhoff stress and

$$\begin{array}{cc}\hfill {\overline{C}}_{ijkl}\left(\overline{\mathit{F}}\right)={\overline{C}}_{ijkl}\left(\overline{\mathit{R}}\overline{\mathit{U}}\right)& =\sum _{m,n=1}^3{\overline{\mathit{R}}}_{im}{\overline{C}}_{mjnl}\left(\overline{\mathit{U}}\right){\overline{\mathit{R}}}_{kn}(i,j,k,l=1,2,3).\hfill \end{array}$$

(42)

for the components of the corresponding stiffness tensor $\mathbb{C}$, see Appendix A. These known facts lead to

Note 2.

In order to collect representative samples of the hyperelastic response functions $\overline{W}$, $\overline{\mathit{P}}$, and $\overline{\mathbb{C}}$, it suffices to evaluate them on samples of the stretch tensor $\overline{\mathit{U}}\sim {\mathbb{R}}^6$ instead of evaluating them on samples of the deformation gradient $\overline{\mathit{F}}\sim {\mathbb{R}}^9$.

This effectively reduces the number of dimensions of $\overline{\mathcal{F}}$ from nine to six. The same dimensionality is attained when considering the response functions with respect to a symmetric measure of strain, e.g., as is done in [30] where the tangent stiffness is efficiently computed using a perturbation technique. However, such measures are unsuitable for reduction by means of the proposed RB method.

The remaining six-dimensional space of s.p.d. tensors is not linear but a convex cone (which does not include the zero element). Moreover, linearly combining elements within this space quickly leads to values of $\overline{J}=det\overline{\mathit{U}}$ describing unphysically large changes in volume. For instance, $\overline{\mathit{U}}=1.3\mathit{I}$ equates to more than 100% volumetric extension, which is well beyond the regime of usual hyperelastic materials that are often nearly incompressible. On the other hand, 100% deviatoric strain is within the range of many standard materials, such as rubber. Hence, in order to describe the space of practically relevant stretch tensors, we propose to apply the DDMS to the macroscopic stretch tensor,

$$\begin{array}{cc}\hfill \overline{\mathit{U}}& ={\overline{J}}^{1/3}\widehat{\overline{\mathit{U}}}.\hfill \end{array}$$

(43)

Let $\widehat{\overline{\mathcal{U}}}$ denote the manifold of unimodular macroscopic stretch tensors $\widehat{\overline{\mathit{U}}}={(\overline{J})}^{-1/3}\overline{\mathit{U}}$. The multiplicative split (43) is the basis for:

Proposition 1.

The set of practically relevant macroscopic stretch tensors $\overline{\mathcal{U}}$ can be sampled via sampling of both the macroscopic determinants,

$$\begin{array}{cc}\hfill \left\{{\overline{J}}^{(m)}{\}}_{m=1}^{N_{\mathrm{dil}}}\right.& \subset {\mathbb{R}}_{+},\hfill \end{array}$$

and the macroscopic unimodular stretch tensors,

$$\begin{array}{cc}\hfill {\left\{{\widehat{\overline{\mathit{U}}}}^{(j)}\right\}}_{j=1}^{N_{\mathrm{dev}}}& \subset \widehat{\overline{\mathcal{U}}},\hfill \end{array}$$

where $N_{\mathrm{dil}}$ and $N_{\mathrm{dev}}$ are the numbers of the samples. The sampling set is determined by the product set

$$\begin{array}{c}\hfill \left\{{\left({\overline{J}}^{(m)}{)}^{1/3}{\widehat{\overline{\mathit{U}}}}^{(j)}\right\}}_{m,j=1}^{m=N_{\mathrm{dil}},j=N_{\mathrm{dev}}}\subset \overline{\mathcal{U}}.\right.\end{array}$$

(44)

The choice of the dilatational samples is relatively straightforward. For many common materials, the expected range of $\overline{J}$ is rather narrow, so a few equisized or adaptive sub-intervals around $\overline{J}=1$ deliver sufficient resolution.

For the space of unimodular s.p.d. matrices representing $\widehat{\overline{\mathit{U}}}\in \widehat{\overline{\mathcal{U}}}$, basic results of Lie group theory can be exploited. We restrict to stating well-known facts that are necessary at this point. For more details, the interested reader is referred to standard text books, such as [31].

Corollary 1.

Let

$$\begin{array}{cc}\hfill \mathrm{symsl}& =\left\{\underline{\underline{Y}}\in {\mathbb{R}}^{3\times 3}|\underline{\underline{Y}}={\underline{\underline{Y}}}^{\mathsf{T}},\mathrm{tr}(\underline{\underline{Y}})=0\right\}\hfill \end{array}$$

be the tangent space and

$$\begin{array}{cc}\hfill {\mathrm{SymSL}}_{+}& =\left\{\underline{\underline{\mathit{U}}}\in {\mathbb{R}}^{3\times 3}|\underline{\underline{\mathit{U}}}={\underline{\underline{\mathit{U}}}}^{\mathsf{T}},det\underline{\underline{\mathit{U}}}=1,{\underline{\mathit{x}}}^{\mathsf{T}}\underline{\underline{\mathit{U}}}\underline{\mathit{x}}>0\forall \underline{\mathit{x}}\in {\mathbb{R}}^3\right\}\hfill \end{array}$$

be the manifold of unimodular s.p.d. matrices. The matrix exponential maps the tangent space bijectively onto the manifold,

$$\begin{array}{c}\hfill \exp :\mathrm{symsl}\to {\mathrm{SymSL}}_{+}.\end{array}$$

The proof of this statement is standard, e.g., by means of the eigenvalue decomposition, and does not necessitate the reference to the abstract setting of Lie groups. In fact, all of the following arguments could be given without the notion of tangent spaces and manifolds. However, this notion is a fundamental concept in nonlinear mechanics. For a very descriptive and comprehensive work on this topic, the reader is referred to [32]. We choose to build upon this concept, as it comes along with vivid interpretations of the function spaces $\overline{\mathcal{U}}$ and $\widehat{\overline{\mathcal{U}}}$.

The set ${\mathrm{SymSL}}_{+}$ is the set of matrix representations of the stretch tensors $\widehat{\overline{\mathit{U}}}\in \widehat{\overline{\mathcal{U}}}$. The tangent space $\mathrm{symsl}$ is the set of Hencky strains, which is linear. Hence, by virtue of the matrix exponential, the sampling of the nonlinear manifold $\widehat{\overline{\mathcal{U}}}$ can be reduced to the sampling of a linear space. Moreover, the restrictions of symmetry and zero trace render the tangent space five-dimensional. This property is, by definition, shared with the manifold ${\mathrm{SymSL}}_{+}$.

The two-dimensional case is now addressed for the sake of visualization. In this setting, the nonlinearity of the manifold and the structure of the sampling set $\overline{\mathcal{U}}$ can be illustrated figuratively. With the subscript (2) denoting two-dimensional quantities, a basis of the tangent space is given by

$$\begin{array}{cccc}\hfill {\underline{\underline{Y}}}_{(2)}^{(1)}& =\sqrt{\frac{1}{2}}\left[\begin{array}{cc}1& 0\\ 0& -1\end{array}\right]\mathrm{and}\hfill & \hfill {\underline{\underline{Y}}}_{(2)}^{(2)}& =\sqrt{\frac{1}{2}}\left[\begin{array}{cc}0& 1\\ 1& 0\end{array}\right].\hfill \end{array}$$

The stretch tensors are obtained through

$$\begin{array}{cc}\hfill {\underline{\underline{\overline{\mathit{U}}}}}_{(2)}& ={\overline{J}}^{\frac{1}{2}}\exp \left(t\left(\alpha {\underline{\underline{Y}}}_{(2)}^{(1)}+\beta {\underline{\underline{Y}}}_{(2)}^{(2)}\right)\right)=\left[\begin{array}{cc}\mathit{a}& \mathit{b}\\ \mathit{b}& d\end{array}\right].\hfill \end{array}$$

A visualization of such samples is depicted in Figure 3. There, for the sake of visual clarity, the determinant $\overline{J}$ is sampled by four equidistant (and rather unrealistic) values between 0.1 and 4. The value $t\in [-2,2]$ is called deviatoric amplitude. A densely uniform sampling $\phi \in [0,2\pi )$ yields the coefficients $\alpha =cos\phi$ and $\beta =sin\phi$.

This emphasizes the important role of the DDMS in the context of sampling, as utilized in (44)—it allows for the identification of a physically meaningful sampling domain that is much smaller than the surrounding space of all admissible stretch tensors. On a side note, the isodet surfaces are perpendicular to the line representing the dilatational stretch tensors.

The proposed sampling procedure for $\overline{\mathcal{U}}$ in three dimensions is given in Algorithm 1. For this purpose, an orthonormal basis ${\underline{\underline{Y}}}^{(1)},\dots ,{\underline{\underline{Y}}}^{(5)}$ of the tangent space $\mathrm{symsl}$ is fixed, cf. Appendix C. The numbers of different kind of samples $N_{\det }$, $N_{\mathrm{dir}}$, and $N_{\mathrm{amp}}$ relate to the quantities $N_{\mathrm{dil}}$ and $N_{\mathrm{dev}}$ introduced in (44) by $N_{\det }=N_{\mathrm{dil}}$ and $N_{\mathrm{dir}}{N}_{\mathrm{amp}}=N_{\mathrm{dev}}$.

Algorithm 1: Sampling of the macroscopic stretch tensor.

Input : ${\overline{J}}_{\min },{\overline{J}}_{\max }$ minimum and maximum determinant with ${\overline{J}}_{\min }<1<{\overline{J}}_{\max }$       $t_{\max }>1$ maximum deviatoric amplitude       $N_{\det }$ number of macroscopic determinants       $N_{\mathrm{dir}}$ number of deviatoric directions       $N_{\mathrm{amp}}$ number of deviatoric amplitudes   Output: $N_{\det }{N}_{\mathrm{dir}}{N}_{\mathrm{amp}}$ samples of $\overline{\mathit{U}}$ 1 Place $N_{\det }$ determinants regularly between the extremal values, $$\begin{array}{c}\hfill {\overline{J}}_{\min }={\overline{J}}^{(1)}<\dots <1<\dots <{\overline{J}}^{\left(N_{\det }\right)}={\overline{J}}_{\max }.\end{array}$$ 2 Generate any approximately uniform distribution of $N_{\mathrm{dir}}$ directions in ${\mathbb{R}}^5$, e.g., as in [12], $$\begin{array}{c}\hfill {\left\{{\underline{N}}^{(n)}\right\}}_{n=1}^{N_{\mathrm{dir}}}\subset \left\{\underline{N}\in {\mathbb{R}}^5:\parallel \underline{N}\parallel =1\right\}.\end{array}$$ 3 Place $N_{\mathrm{amp}}$ deviatoric amplitudes regularly between 0 and the expected maximum value, $$\begin{array}{c}\hfill 0\le t^{(1)}<\dots <t^{\left(N_{\mathrm{amp}}\right)}=t_{\max }.\end{array}$$ 4 Return the set of samples of $\overline{\mathit{U}}$: $$\begin{array}{c}\hfill \left\{{\left({\overline{J}}^{(m)}{)}^{1/3}\exp \left(t^{(p)}\left[\sum _{k=1}^5{\left({\underline{N}}^{(n)}\right)}_k{\underline{\underline{Y}}}^{(k)}\right]\right)\right\}}_{m,n,p=1}^{m=N_{\det },n=N_{\mathrm{dir}},p=N_{\mathrm{amp}}}\subset \overline{\mathcal{U}}\right.\end{array}$$ (45)

The order of Steps 1 to 3 is interchangeable. Details on these parts are now presented:

Step 1.

Uniform seeding of the determinants is actually not required, but any pattern implying the sampling determinants ${\{{\overline{J}}^{(k)}\}}_{k=1}^{N_{\det }}$ to be dense in $[{\overline{J}}_{\min },{\overline{J}}_{\max }]$ as $N_{\det }\to \infty$ works without loss of generality. In this way, the dilatational response may be resolved adaptively.

Step 2.

The generation of uniform point distributions on spheres is a research topic on its own, see [33] for an overview. The method described in [12] is based on energy minimization, which is also used in the present work. Some point sets of various sizes are included in the example program [18]. More detailed investigations on this topic and an open-source code of a point generation program are part of another work, [34]. Alternatively, Equal Area Points [35] may be used as a rough but quickly computable approximation of such point sets.

Step 3.

As in Step 1, the uniform placement of the deviatoric amplitudes, $t^{(p)}$, may be substituted by adaptive alternatives. In [12], we have suggested to use an exponential distance function.

The result of Steps 2 and 3, i.e., the sampling of the tangent space, is exemplified in Figure 4 for the two-dimensional case ($\mathit{u}\in {\mathbb{R}}^2$) and for $N_{\mathrm{dir}}=5$ and $N_{\mathrm{amp}}=3$. There, an adaptive spacing of the deviatoric amplitudes is applied. This might be beneficial for capturing strongly changing material behavior near the relaxed state.

In general, the vector $\underline{N}\in {\mathbb{R}}^5$, $\parallel \underline{N}\parallel =1$, corresponding to a macroscopic Hencky strain characterizes the direction of the applied stretch $\overline{\mathit{U}}$, which we also coin the load case.

3.3. Application of the Stretch Tensor Trained Reduced Basis Model

Since the RB model is trained on only the rotationally invariant part of $\overline{\mathit{F}}$ but should be applied to general deformation gradients, the transformation rules (41) and (42) are incorporated during the evaluation of the surrogate model. Details on the procedure are given in Algorithm 2.

Algorithm 2: Online phase of the stretch tensor trained Reduced Basis method

Input : $\overline{\mathit{F}}$ macroscopic deformation gradient   Output: ${\overline{\mathit{P}}}^{\mathrm{RB}}(\overline{\mathit{F}})$, ${\overline{\mathbb{C}}}^{\mathrm{RB}}(\overline{\mathit{F}})$ effective material response 1  Compute polar decomposition $\overline{\mathit{F}}=\overline{\mathit{R}}\overline{\mathit{U}}$. 2  Compute approximations of effective stress ${\overline{\mathit{P}}}^{\mathrm{RB}}(\overline{\mathit{U}})$ and effective stiffness ${\overline{\mathbb{C}}}^{\mathrm{RB}}(\overline{\mathit{U}})$. 3  Transform effective responses to correct frame ${\overline{\mathit{P}}}^{\mathrm{RB}}(\overline{\mathit{F}})$, ${\overline{\mathbb{C}}}^{\mathrm{RB}}(\overline{\mathit{F}})$, using $\overline{\mathit{R}}$, cf. (41), (42).

4. Numerical Examples

4.1. Reduced Basis for a Fibrous Microstructure

The applicability of the proposed RB method in combination with the sampling scheme is now numerically studied for a fibrous microstructure roughly resembling polymers with woven reinforcements. The goal is to prove the efficiency of the F-RB scheme in principle and under “worst-case” conditions, the latter meaning that the phase contrast is chosen to be rather large. Yet, at this stage, it is explicitly not aspired to provide fully realistic examples. These would require investigations on the proper size of the microstructure and should employ dissipative material laws, both of which are not within the scope of this work.

A cubical microstructure with two fibrous inclusions is considered, see Figure 5a and cf. [36] for a related example. The inclusions are periodic and occupy approximately 0.7% of the volume. The mesh contains 35,516 nodes in 25,633 quadratic tetrahedron elements (C3D10).

For the matrix, the material constants are chosen to be $K_{\mathrm{m}}=400$ MPa and $G_{\mathrm{m}}=0.4$ MPa, resembling rubber-like material properties. For the fibers, the values are $K_{\mathrm{f}}=800$ MPa and $G_{\mathrm{f}}=240$ MPa. The latter parameters approximate the behavior of stiffer polymers, such as polyethylene. The phase contrast is 600 with respect to the shear moduli, and the Poisson ratios are ${\nu }_{\mathrm{m}}=0.4995$ and ${\nu }_{\mathrm{f}}=0.3636$.

The training boundary conditions are defined with $N_{\mathrm{dir}}=128$ deviatoric directions, ${\underline{N}}^{(n)}$, and $N_{\mathrm{amp}}=10$ deviatoric amplitudes, $t^{(p)}$, which are regularly distributed in the interval $[0.05,0.5]$. In order to also consider response to volumetric extension in the training data, an additional set of $N_{\det }$ = 10 boundary conditions of the form $({\overline{J}}^{(m)}{)}^{1/3}\mathit{I}$ is included in the training set, with the determinant ${\overline{J}}^{(m)}$ being linearly increased from 1 to 1.02.

The validation load cases are 640 mixed dilatational-deviatoric boundary conditions. Along $N_{\mathrm{dir}}=64$ new deviatoric directions, both the deviatoric amplitudes ($t^{(p)}=0.05,\dots ,0.5$) and the dilatational amplitudes (${\overline{J}}^{(m)}=0.9995,\dots ,0.995$) are applied in 10 equidistant increments.

The results for various values N of the RB-size are compared with the results of FE simulations with the same boundary conditions. To this end, the error measures

$$\begin{array}{ccccc}\hfill {\mathrm{err}}_{\mathrm{W}}& =\frac{\parallel {\overline{W}}^{\mathrm{RB}}-{\overline{W}}^{\mathrm{FEM}}\parallel }{\parallel {\overline{W}}^{\mathrm{FEM}}\parallel }\hfill & \mathrm{and}& \hfill {\mathrm{err}}_{\mathrm{P}}& =\frac{\parallel {\overline{\mathit{P}}}^{\mathrm{RB}}-{\overline{\mathit{P}}}^{\mathrm{FEM}}\parallel }{\parallel {\overline{\mathit{P}}}^{\mathrm{FEM}}\parallel }\hfill \end{array}$$

(46)

are employed. Figure 6 and Figure 7 visualize the results.

The distribution of the energy error, ${\mathrm{err}}_{\mathrm{W}}$, improves monotonically as the RB is enriched from $N=8$ to $N=128$ elements. This enrichment corresponds to the inclusion of additional subtrahends in the computation of ${\overline{\mathbb{C}}}^{\mathrm{RB}}$, improving the spectral over-estimation by the volume average of the stiffness, cf. (17). It is also noteworthy that the error tends to be higher for larger magnitudes of the applied kinematic boundary condition, although that is not always the case.

In contrast to this, the stress error ${\mathrm{err}}_{\mathrm{P}}$ distribution monotonically improves only up to a certain threshold value of the number of basis elements, which is $N=64$ in this example. For the greater bases with $N=96$ and $N=128$, the quality of the results deteriorates as far as the stress error is concerned. This is most likely due to an excessively oscillatory nature of the higher order modes—at some critical level $1\ll i=N_{\mathrm{crit}}<N$, the POD constructs eigenvectors ${\underline{E}}^{(i)}$ with the $L^2({\Omega }_0)$-norm $\sqrt{{\lambda }_i}\ll 1$. Therefore, the POD would construct basis vectors out of numerical fluctuations, which would be unphysical. While the enrichment of the optimization space with unphysical information cannot increase the minimum energy error ${\mathrm{err}}_{\mathrm{W}}$, it might lead to fluctuations in the displacement field that have significant impact on the overall stress response. This is especially the case for numerical fluctuations within the stiff inclusion phase where low overall displacement errors still could lead to notable impact on the effective stress.

Nonetheless, all observed errors are less than 20% and stay below 3% for the optimal sampled size $N=64$. For half the basis size, $N=32$, the errors max at approximately 5%, which is still acceptable considering the uncertainties involved in realistic two-scale simulations. Note that the maximum errors strongly depend on the maximum load amplitude, which is chosen to be very large in this example (50% deviatoric strain and 0.5% compression).

The runtimes of the RB model for different sizes N are depicted in Figure 8. A nearly linear growth of the runtime with respect to the basis size can be asserted. It is noteworthy that the online time of the RB method is strongly dominated by the assembly of the Jacobian $\underline{\underline{D}}$. Therefore, a Quasi-Newton implementation was chosen, resulting in only two assemblies per load increment.

Speed-ups become impressive when very large load increments are considered. In all examples observed thus far, the RB converges to the final load amplitude in a single increment, requiring 7–13 Quasi–Newton iterations, with only 2–4 assemblies of the Jacobi matrix $\underline{\underline{D}}$ and a runtime of 10–50 seconds. This is in strong contrast to the FEM which is very sensitive to large load increments as they come along with a high probability of a violation of the condition $\det (\mathit{F})>0$. By means of standard implementation, such occurrences are usually treated with cutbacks of the load increment, which is detrimental to the runtime of the FEM.

No rigorous speed-up analysis is intended at this point. Both the codes of the FEM and of the RB method are fairly efficient in-house C/C++ developments and perform exact line searches. While the FEM has not yet been equipped with a Quasi–Newton procedure, the linear solver makes use of parallelization. This is in contrary to the current implementation of the RB method. Depending on the microstructure (especially the geometry, material nonlinearities, and phase contrast), the loading conditions, and the size N of the RB, speed-up factors are in the order of 5–100.

4.2. Reduced Basis for a Stiffening Microstructure

The second example takes the “worst-case” approach further by considering a noncubical microstructure with even higher phase contrast and significant topological nonlinearity. To this end, a cuboid microstructure occupying the domain $[-0.5,0.5]\times [-0.3,0.3]\times [-0.05,0.05]\subset {\mathbb{R}}^3$ and containing a hash-like inclusion is investigated, see Figure 9a. The mesh is periodic and contains 33,923 nodes in 21,726 quadratic tetrahedron elements (C3D10). The reinforcement makes up approximately 13.3% of the volume. Due to this large volume fraction, a pronounced geometry-induced nonlinearity of the effective response is expected under uniaxial loading conditions along the x-axis. As it is elongated, the hashlike part is straightened and thus increasingly aligned with the external loading, see Figure 9b. Such effects might be desirable when designing microstructures for functional materials.

The material parameters are $K_{\mathrm{m}}=19.867$ MPa, $G_{\mathrm{m}}=0.4$ MPa, $K_{\mathrm{f}}=\mathrm{19,867}$ MPa, and $G_{\mathrm{f}}=400$ MPa, implying a Poisson ratio of $0.49$ in both materials and a phase contrast of 1000.

The training boundary conditions are the deviatoric ones of the set considered in Section 4.1, i.e., $N_{\mathrm{dir}}=128$ deviatoric directions and $N_{\mathrm{amp}}=10$ regularly spaced deviatoric amplitudes from the interval $[0.05,0.5]$. No dilatational training cases are considered, i.e., only points from a five-dimensional submanifold of the space $\overline{\mathcal{U}}$ are sampled.

Uniaxial tension boundary conditions are applied for the validation. More precisely, the stretch component ${\overline{\mathit{U}}}_{\mathrm{xx}}$ is increased from $1.0$ to $1.5$ in 10 increments of equal size. The other components are chosen such that all but the $\mathrm{xx}$-component of the effective stress $\overline{\mathit{P}}$ vanish.

Figure 10 depicts the results for different sizes N of the RB. The influence of the stiffening effect on the stress curve is emphasized by the black dashed line corresponding to a similar microstructure with a straight, cuboid inclusion that leads to the same final stress value under these boundary conditions, see Figure 9c.

In this example, the geometric stiffening effect is captured by the RB with high accuracy, with as few as $N=24$ basis elements. For moderate stretches, even an RB size of $N=16$ suffices. These results are somewhat more impressive when noticing that the applied boundary condition contains more than 1.2% volumetric compression, i.e. the validation loading actually lies outside the submanifold covered during training.

In order to examine the action of the cutoff function ϕ, the following two indicators are introduced:

$$\begin{array}{cc}\hfill c_{\mathrm{qp}}& =\#\mathrm{quadrature}\mathrm{points}\mathrm{with}\left(\varphi (J)<1\right),\hfill \end{array}$$

(47)

$$\begin{array}{cc}\hfill V_{\mathrm{excl}}& =\left(|{\Omega }_0|-\sum _{p=1}^{N_{\mathrm{qp}}}\varphi (J_p){\mathit{w}}_p\right)/|{\Omega }_0|.\hfill \end{array}$$

(48)

The first quantity, $c_{\mathrm{qp}}$, counts the number of quadrature points at which the cutoff function has an influence. The second one, $V_{\mathrm{excl}}$, measures the relative excluded volume, interpreting the value of ϕ as a scaling of the corresponding quadrature weight. The values of these indicators are depicted in Figure 11.

Most notably, the cutoff function does not have any effect before the eighth load increment in this example. Only for large load amplitudes does this numerical stability tweak become necessary. Even then, the number of points at which it has an influence is small, considering the total number of quadrature points, $N_{\mathrm{qp}}=\mathrm{86,904}$. This example is representative for all conducted numerical experiments.

5. Discussion

5.1. Discussion of the Reduced Basis Method

5.1.1. Relation of the RB Homogenization to Analytical Estimates

Zero coefficients, $\underline{\xi }=\underline{0}$, correspond to the Taylor homogenization, i.e., to the nonlinear counterpart of the Voigt estimate [37], which provides an upper bound for the material response, cf. (17). Starting with the initial guess ${\underline{\xi }}^{(0)}=\underline{0}$, the evolution of the coefficients corresponds to a (possibly not monotonous) relaxation of this overly stiff response into a more natural state. In view of improved computational efficiency, a nonzero initial guess ${\underline{\xi }}^{(0)}$ combined with an exact line search has proven reasonable and easy to realize. For instance, such a guess might stem from previous load steps or an interpolation/extrapolation of available coefficient data.

5.1.2. Reconstruction of Displacement Fields

Given an RB approximation of the deformation gradient, ${\mathit{F}}^{\mathrm{RB}}$, one can reconstruct the corresponding displacement field uniquely up to rigid body motion. This is possible due to the linear dependence of the deformation gradient fluctuations on the displacement fluctuations. Recall the definition of the RB in (32),

$$\begin{array}{cc}\hfill {\mathit{B}}^{(i)}& =\sum _{j=1}^{N_{\mathrm{s}}}\frac{1}{\sqrt{{\lambda }_i}}{\left({\underline{E}}^{(i)}\right)}_j{\tilde{\mathit{F}}}^{(j)}(i=1,\dots ,N).\hfill \end{array}$$

The corresponding displacement fluctuations are

$$\begin{array}{cc}\hfill {\tilde{\mathit{u}}}_{\mathrm{B}}^{(i)}& =\sum _{j=1}^{N_s}\frac{1}{\sqrt{{\lambda }_i}}{\left({\underline{E}}^{(i)}\right)}_j{\tilde{\mathit{u}}}^{(j)}(i=1,\dots ,N).\hfill \end{array}$$

(49)

The displacement fluctuation fields ${\tilde{\mathit{u}}}^{(j)}$ are defined by ${\tilde{\mathit{u}}}^{(j)}(\mathit{X})={\mathit{u}}^{(j)}(\mathit{X})-{\overline{\mathit{H}}}^{(j)}\mathit{X}$, where the displacement fields ${\mathit{u}}^{(j)}$ are the solutions computed during training, and ${\overline{\mathit{H}}}^{(j)}={\overline{\mathit{U}}}^{(j)}-\mathit{I}$. Thus, a displacement field compatible to the RB result ${\mathit{F}}^{\mathrm{RB}}(\mathit{X};\overline{\mathit{F}})$ is given by

$$\begin{array}{cc}\hfill {\mathit{u}}^{\mathrm{RB}}(\mathit{X};\overline{\mathit{F}})& =\overline{\mathit{H}}\mathit{X}+\sum _{i=1}^N{\xi }_i^{*}(\overline{\mathit{F}}){\tilde{\mathit{u}}}_{\mathrm{B}}^{(i)}(\mathit{X}).\hfill \end{array}$$

(50)

The missing term $\overline{\mathit{u}}(\overline{\mathit{X}})$, cf. (10), cannot be reconstructed.

5.1.3. Relation to Classical Displacement-Based POD Methods

In a certain sense, the entries of the correlation matrix used in the offline phase, cf. Section 2.2, are “weighted” scalar products of the displacement fluctuation fields ${\tilde{\mathit{w}}}^{(i)}$ within the Sobolev space ${\mathit{H}}_0^1({\Omega }_0)$. “Weighted” is to be understood in that the zeroth order derivative is multiplied by zero. Classical displacement-based POD methods compute correlations of the fluctuations ${\tilde{\mathit{w}}}^{(i)}$ within the Lebesgue space $L^2({\Omega }_0)$. The change to ${\mathit{H}}_0^1({\Omega }_0)$-like scalar products is physically motivated by the fact that the local energy $W=W(\mathit{F})$ explicitly depends on the gradient of the displacement, $\mathit{F}=\mathit{u}\otimes {\mathsf{\nabla }}_{\mathrm{X}}+\mathit{I}$, but not on the displacement, u.

5.1.4. Advantages Compared to General Displacement-Based Schemes

The proposed method is advantageous compared to both displacement-based POD methods and the classical FEM for the following reasons:

• No gradients need to be computed from displacement fields, which displacement-based schemes always require prior to the evaluation of the material law.
• The residual $\underline{r}$ and the Jacobian $\underline{\underline{D}}$ are algorithmically sleek and trivial to implement.
• The absence of element formulations in the assembly of the reduced residual $\underline{r}$ and of the Jacobi matrix $\underline{\underline{D}}$ contributes to both the simplicity and the efficiency of the method—no incidence matrices occur, allowing for linear memory access. Moreover, the algebraic operations associated with reference element formulations are bypassed. This is also in favor of parallel computations. Such an implementation is still outstanding for the problem at hand, but has been conducted for related problems in the small strain setting in [38].

Although the storage of the basis ${\{{\mathit{B}}^{(i)}\}}_{i=1}^N$ requires $9N_{\mathrm{qp}}N$ double precision values, the basis is compact enough to be completely loaded into random access memory of standard computers. For instance, the bases considered in Section 4 occupy only ~200 Mb of memory for $N=32$.

We now briefly address the algorithmic complexity associated with the proposed F-RB method and with the u-RB method that was employed in previous works, such as [9] and [8]. To this end, the fully discretized versions of the residual $\underline{r}$ and of the Jacobian $\underline{\underline{D}}$ as well as discrete quantities associated with the u-RB method are introduced in the following listing.

• $\underline{\mathit{P}}({\mathit{X}}_p)\in {\mathbb{R}}^9$: Nine values of the stress $\mathit{P}({\mathit{F}}_{\xi }({\mathit{X}}_p))$ at the quadrature point ${\mathit{X}}_p\in {\Omega }_0$
• $\underline{\underline{C}}({\mathit{X}}_p)\in {\mathbb{R}}^{9\times 9}$: Symmetric stiffness tensor
• $\underline{\underline{\mathit{B}}}({\mathit{X}}_p)\in {\mathbb{R}}^{9\times N}$: F-RB matrix containing the nine values of each basis elements ${\mathit{B}}^{(i)}$ as columns
• ${\mathit{w}}_p$: The quadrature weight at ${\mathit{X}}_p$
• $N_{\mathrm{DOF}}$: Three times the number of nodes, $\sim N_{\mathrm{qp}}$
• ${\underline{r}}^{\mathrm{FE}}\in {\mathbb{R}}^{N_{\mathrm{DOF}}}$: Global FE residual vector
• ${\underline{\underline{\mathit{B}}}}^{\mathrm{u}}\in {\mathbb{R}}^{N_{\mathrm{DOF}}\times N}$: u-RB matrix of which the columns contain the nodal displacement values
• ${\underline{\underline{K}}}^{\mathrm{FE}}\in {\mathbb{R}}^{N_{\mathrm{DOF}}\times N_{\mathrm{DOF}}}$: Global FE stiffness matrix

Table 1. Algorithmic complexities of the well-established u-based RB method and the novel F-based RB method. The assembly of the FE residual and of the FE stiffness depend on $N_{\mathrm{qp}}$.

RB Method	Quantity	Complexity
F-based	${\displaystyle \underline{r}=\sum _{p=1}^{N_{\mathrm{qp}}}\underline{\underline{\mathit{B}}}{({\mathit{X}}_p)}^{\mathsf{T}}\underline{\mathit{P}}({\mathit{X}}_p){\mathit{w}}_p}$	$\mathcal{O}(9N{N}_{\mathrm{qp}})$
	${\displaystyle \underline{\underline{D}}=\sum _{p=1}^{N_{\mathrm{qp}}}\underline{\underline{\mathit{B}}}{({\mathit{X}}_p)}^{\mathsf{T}}\underline{\underline{C}}({\mathit{X}}_p)\underline{\underline{\mathit{B}}}({\mathit{X}}_p){\mathit{w}}_p}$	$\mathcal{O}([9^{2}N+9N^2]N_{\mathrm{qp}})$
u-based	${\displaystyle \underline{r}={\left({\underline{\underline{\mathit{B}}}}^{\mathrm{u}}\right)}^{\mathsf{T}}{\underline{r}}^{\mathrm{FE}}}$	$\mathcal{O}(N{N}_{\mathrm{DOF}})$ + assembly of ${\underline{r}}^{\mathrm{FE}}$
	${\displaystyle \underline{\underline{D}}={\left({\underline{\underline{\mathit{B}}}}^{\mathrm{u}}\right)}^{\mathsf{T}}{\underline{\underline{K}}}^{\mathrm{FE}}{\underline{\underline{\mathit{B}}}}^{\mathrm{u}}}$	$\mathcal{O}(N{N}_{\mathrm{DOF}}^2+N^2{N}_{\mathrm{DOF}})$ + assembly of ${\underline{\underline{K}}}^{\mathrm{FE}}$

compares the algorithmic complexity of the presented F-RB method with that of standard u-RB schemes. First of all, both methods share a quadratic dependence of their Jacobi matrices on the number of modes, N. Therefore, the assembly of the Jacobian is usually the most costly operation. Secondly, both approaches’ complexities suffer a linear dependency on the number of quadrature points, $N_{\mathrm{qp}}$. In the displacement-based approach, this is included in the assembly of the residual and of the stiffness, which relate to the factor 9 and 9², respectively. Thirdly, the novel F-RB scheme spares the computational overhead associated with FE formulations ${\underline{r}}^{\mathrm{FE}}$ and ${\underline{\underline{K}}}^{\mathrm{FE}}$. More details on this matter are currently being investigated.

5.1.5. Outlook

Future research should aim at an application of the introduced Reduced Basis method within realistic two-scale simulations, in analog to [12,38,39,40]. Hyperreduction methods, cf. [41], might give rise to additional speed-ups in the online phase. Further, modifications of the cutoff function, ϕ, should be investigated—a function with compact support might be more appropriate. The construction of the RB from large sets of snapshots is computationally intense, as much data needs to be processed. In the above examples, the POD consumed multiple hours of time. Hierarchical approximations, such as [42], might mitigate the effects by enabling parallel computations. Overall, the long-term perspective is to extend this RB framework efficiently to the context of dissipative materials.

5.2. Discussion of the Sampling Strategy

The proposed sampling strategy is meant to serve as a template. As exemplified in Section 4.1, the samplings can be modified and still lead to a coverage of the set of macroscopic boundary conditions that is sufficient for the problem at hand. The example of Section 4.2 took this idea further and showed that it might not even be necessary to sample the macroscopic determinant. Hence, the sampling can sometimes be reduced to the five-dimensional subspace of isochoric macroscopic stretch tensors.

In any case, the exact choice of both the inputs to Algorithm 1 and the distributions of the deviatoric amplitudes and the macroscopic determinants remains to be based on knowledge and sophisticated guesses, at least at the current state of the art. Further research on this matter might lead to a refined alternative to Algorithm 1, possibly involving the evaluation of error estimators.