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, in which 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 extend. Therefore, it is worthwile 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 (Rendek and Lion, 2010) and amorphous glassy polymers (Nguyen et al., 2016), 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. 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 Geers and Yvonnet (2016) and Saeb et al. (2016), summarize the progress. Here, attention is restrained to few methods most similar or relevant to the present work.
The FE2 method Feyel (1999) 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 FE2, 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. Such methods include the R3M (Yvonnet and He, 2007) and can be further enhanced by the use of, e.g., the EIM, as in Radermacher and Reese (2016). Numerical comparisons of various schemes were conducted in Radermacher and Reese (2013) and Soldner et al. (2017)
. 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.
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 is made in the field of efficient integration schemes, see for instance An et al. (2009) and Hernández et al. (2017). 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 Temizer and Zohdi (2007)interpolates the coefficients of an assumed macroscopic material model. Another example is the NEXP method Yvonnet et al. (2013), 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 Fritzen and Kunc (2018), 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 the 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 work Fritzen and Kunc (2018) to the finite strain regime. It aims at reducing the computational complexity of the problem associated with the microstructure 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 demonstrator containing a minimum working example of the RB model with 50 lines of MATLAB/Octave code, (Kunc, 2019).
As for the setup phase, the snapshot data is created by means of an efficient sampling procedure. To this end, a linear parameter space is identified, within which the sampling sites are placed based upon physical interpretation. This allows for controllability of the resolution of certain key characteristics of the effective material response while keeping the total number of samples within bounds.
The set of real numbers and the subset of positive numbers greater than zero are denoted by and
, respectively. Matrices are marked by two underlines and vectors by one underline, e.g., . Vectors are assumed to be columns and the dot product of two vectors of the same size is understood as the Euclidean scalar product, . First order and second order tensors in coordinate free description are denoted by bold letters, e.g. , . No conclusion can be drawn on the order of a tensor based on its capitalization. Here, the underlying space is always the Euclidean space with its standard basis. First order and second order tensors can be represented as vectors and matrices, e.g. and 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 , e.g. and . Components of tensors of order are with respect to the Euclidean tensorial basis , e.g. , for second order tensors , and , for , . The following contractions are defined:
Let be the domain occupied by a physical body undergoing elastic deformations, and let be its initial configuration. Then and describe the coordinates of material points within the current configuration and within the reference state respectively. Their difference is the displacement , see Figure 1. The gradient of a vector field is defined as a right gradient and denoted by . 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 and . The set of square-integrable Lebesgue functions on the reference domain is tagged .
The displacement gradient and the deformation gradient are related through , where is the second order identity tensor in three dimensions. The determinant measures the relative volumetric change due to the present deformation.
Unimodular quantities, i.e. second order tensors with determinant one, may be emphasized by a hat, e.g. . This multiplicative decomposition is sometimes attributed to Flory Flory (1961) 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. , , while symbols without overline correspond to their microscopic counterpart, e.g. , . Equivalently, macroscopic quantities are called global and microscopic ones are called local. The volume average of a general local field
is essential to the theory. The dependence of a microscopic quantity on both the microscopic coordinates and a macroscopic quantity is denoted by . In such a case, the components of the macroscopic quantity are called parameters of the microscopic function . The application of the volume averaging operator is abbreviated by . The case of a concatenated function is analogous, i.e. , regardless of the tensorial order of the image of the function .
1.5 Material models
In this work hyperelastic materials are investigated. They are characterized by stored energy density functions . The first Piola-Kirchhoff stress
|and the corresponding fourth-order stiffness tensor|
characterize the material response.
Henceforth, for reasons of readability, the stored energy density function 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 that model deformation plasticity (e.g., Fritzen and Kunc, 2018; Yvonnet et al., 2013; Bilger et al., 2005). Although these models are only valid for purely proportional loading conditions, they provide means to simulate highly nonlinear material behavior in certain scenarios comparably easy 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
is used, with the bulk modulus and the shear modulus. The volumetric part of the energy is taken from Doll and Schweizerhof (1999). 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. (Simo, 1988).
1.6 Problem setting of first order homogenization
Assuming stationarity and separability of scales, the following coupled, deformation-driven problems can be derived by means of asymptotic expansion of the displacement and subsequent first order approximation. This procedure is carried out in Pruchnicki (1998) 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
where denote bulk forces, and balance of angular momentum
along with well-posed boundary conditions constitute the macroscopic boundary value problem. This system of equations is closed by means of the hyperelasticity law, cf. (2),
Note that, in general, is a priori not available and is thus a purely formal relation. For reasons of readability, the dependence of any quantity on the macroscopic material coordinate is usually spared, e.g. , , or .
1.6.2 Microscopic problem
The microscopic boundary value problem is given by the balance equations
and suitable boundary conditions, e.g. as discussed in Miehe (2003). In this work, periodic displacement fluctuation boundary conditions are employed. The microscopic displacements take the form
Therein, the macroscopic displacement is independent of microscopic quantities. The second term, , corresponds to a homogeneous deformation of the microstructure. The third term, , is a displacement fluctuation with the zero mean property . Hence, the deformation gradient reads
where the fluctuational part , has the zero mean property
Thus, the volume average of the local deformation gradient equals its macroscopic counterpart,
This motivates the identification of as the boundary condition to the microscopic problem (8). As for the material response, the following relationships can be deduced:
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,
2 Reduced Basis homogenization for hyperelasticity
The Reduced Basis (RB) scheme is based on a direct approximation of the microscopic deformation gradient from Equation (11) without the need to explicitly have the corresponding displacement at hand. The initial approximation is given by
The parameters and are the boundary condition to (8) and the reduced coefficient vector, respectively. Note that the macroscopic coordinates is not assumed to be an explicit part of the parameter set, i.e. the same approximation is made throughout the macrostructure. The set is linearly independent within the space and is called RB of . In a later context, it will also be referred to as the set of ansatz functions. In order to enforce the relationship
regardless of the reduced coefficients , the basis functions are asserted to have the fluctuation property, i.e. for
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 and of the reduced coefficients ,
By the principle of minimum energy, the optimal coefficients
are sought-after. The unconstrained and unique solvability of this task is assumed for the moment and will be addressed in Section2.4. The solution of (22) defines the RB approximation of the deformation gradient
The microscopic energy, stress and stiffness within the microstructure are then approximated by
respectively. The resulting effective energy is readily given by
Also, the effective responses and 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 which 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 ,
The solution stress field is -orthogonal to the RB ansatz functions .
A viable choice for the solution of the minimization problem (22) is the Newton-Raphson scheme, which necessitates the Jacobian with the components
Then, the th iteration to the solution reads
The initial guess can be zero or a more sophisticated alternative.
In total, the problem of determining both the local field and the homogenized material properties depends only on degrees of freedom, namely on the coefficients . Usually, the number of DOF is in the range of 10 to 150, which compares to the full order model’s complexity that can easily reach more than 105 DOF.
Despite this impressive reduction of the number of DOF, the computational costs associated with the assembly of the residual and of the Jacobian 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 is considered. For more information on these topics, the reader is referred to the authors’ previous work Fritzen and Kunc (2018) or standard literature, such as (Castañeda and Suquet, 1998, part II.C).
2.2 Identification of the Reduced Basis
The basis is computed by means of a classical snapshot POD, see (Sirovich, 1987). 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. Instead, the POD is performed on fluctuations of the deformation gradient.
During the snapshot POD, first, snapshots are created by means of high-fidelity solutions to the nonlinear microscopic problem (8) with different snapshot boundary conditions , , which are also called training boundary conditions. Each of these boundary conditions leads to a solution field
. Typically, the Finite Element Method (FEM) or solvers making use of the Fast Fourier Transform(e.g. Kabel et al., 2014) 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 . 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,
Each of these fluctuation fields represents 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 scalar products of the snapshots, (). The remaining procedure is standard, see for instance (Sirovich, 1987) or (Quarteroni et al., 2016): the eigenvalues
, corresponding to the eigenvectors, are sorted in descending order and truncated by only considering the first values, . The decision on a particular threshold index is based on consideration of the Schmidt-Eckhard-Young-Mirsky theorem. Finally, the RB is constructed as
where the factor accounts for normalization of the basis elements. We conclude that the RB is a collection of orthonormal -like quantities.
2.3 Mathematical motivation of the Reduced Basis model
Next, the obtained deformation gradient field and the associated stress field are shown to weakly solve the original problem (8), .
Let be an admissible test function, i.e. a once weakly differentiable periodic displacement fluctuation field, and let denote its gradient. Then the well-known weak form of (8) is equivalent to the principle of virtual work,
The residual from (28) coincides with the left-hand side of the weak form, if the test function is chosen suitably. In fact, for each basis element there is a corresponding displacement fluctuation field that is defined uniquely up to rigid body motion, hence the equivalence of (28) and (33).
It follows that the Reduced Basis scheme is a Galerkin procedure, taking the displacement fields corresponding to the RB elements as both ansatz and test functions. Hence, the RB model is equivalent to the FEM, but with basis functions that are globally supported in . In other words, the basis functions of the RB method span subspace of dimension 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 Ball (1976), and recent developments in this matter, e.g. Schneider (2017).
A constraint to the optimization problem is the physical condition
which means that no singular () or self-penetrating () deformations shall occur. This reduces the set of admissible coefficients to a subset of that is non-trivial 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
Here, , , and 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 might possibly have some positive but overly small values of the determinant, , 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 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:
Therein, the almost smooth cutoff function is empirically defined by
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 vectorsfor which effective quantities can be computed, albeit approximately, to the whole space .
The significance of local fields varies with the value of the cutoff function. When attains values less than one, information is considered accordingly less reliable. In this sense, microscopic information is filtered based on a trust region for defined by can be seen as a reliability indicator.
3.1 General considerations
The proposed sampling strategy builds on the previous work Fritzen and Kunc (2018), 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, , of the macroscopic deformation gradient, , and the nonlinearity of the volumetric part of the deformation, , with respect to the local displacements, . For the setup of the Reduced Basis model as described in Section 2.2, the space of macroscopic deformation gradients
needs to be sampled, i.e. the discrete sampling set is to be defined. Two contradictory requirements need to be satisfied when constructing :
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 wholistic 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 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 is a subset of a nine-dimensional space () restricted by the inequality
Therefore, a regular grid in the components of 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 solutiotns of the full order model. 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, , where is a rotation and is the symmetric positive definite (s.p.d.) stretch tensor. Material objectivity implies the energy function to be independent of the frame of reference,
|and the transformation behavior|
|for the first Piola-Kirchhoff stress and|
for the components of the corresponding stiffness tensor , see Appendix A. These known facts lead to the following
In order to collect representative samples of the hyperelastic response functions , , and , it suffices to evaluate them on samples of the stretch tensor instead of evaluating them on samples of the deformation gradient .
This effectively reduces the number of dimensions of 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 Miehe (1996) 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 describing unphysically large changes in volume. For instance, , with the second-order identity tensor , 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,
Let denote the manifold of unimodular macroscopic stretch tensors . The multiplicative split (43) is the basis for The set of practically relevant macroscopic stretch tensors can be sampled via sampling of both the macroscopic determinants,
|and the macroscopic unimodular stretch tensors,|
The sampling set is determined by the product set
The choice of the dilatational samples is relatively straight-forward. For many common materials, the expected range of is rather narrow, so a few equisized or adaptive sub-intervals around deliver sufficient resolution.
For the space of unimodular s.p.d. matrices representing , 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 Faraut (2008).
|be the tangent space and|
be the manifold of unimodular s.p.d. matrices. The matrix exponential maps the tangent space bijectively onto the manifold,
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 refered to Neff et al. (2016). We choose to build upon this concept, as it comes along with vivid interpretations of the function spaces and .
The set is the set of matrix representations of . Notably, the tangent space is linear. Hence, by virtue of the matrix exponential, the sampling of the nonlinear manifild 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 .
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 can be illustrated figuratively. With the subscript (2) denoting two-dimensional quantities, a basis of the tangent space is given by
The stretch tensors are obtained through
A visualization of such samples is depicted in Figure 3. There, for the sake of visual clarity, the determinant is sampled by four equidistant (and rather unrealistic) values between 0.1 and 4. The value is called deviatoric amplitude. A densely uniform sampling yields the coefficients and .
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 in three dimensions is given in Algorithm 1. For this purpose, an orthonormal basis of the tangent space is fixed, cf. Appendix C. The numbers of different kind of samples , , and relate to the quantities and introduced in (44) by and .
Generate any approximately uniform distribution ofdirections in , e.g. as in Fritzen and Kunc (2018),
The order of Steps 1 to 3 is interchangeable. Details on these parts are now presented:
Uniform seeding of the determinants is actually not required, but any pattern implying the sampling determinants to be dense in as