1 Introduction
Homogenization in engineering
In the last decade, the quality of micro xray computed tomography (CT) images has steadily improved. Nowadays, standard CT devices have a spatial resolution below one m. They produce 3D images of up to 4096^{3} voxels. This permits a detailed view of the microstructure’s geometry of composite materials all the way down to the limits of continuum mechanical theories and the necessity for discrete particle methods. The geometrical information itself is sufficient to detect defects in components by measuring, e.g., the pore space and the shape of the pores or by detecting the presence of microcracks. In order to understand the effect of such microscopic features, one has to solve PDEs on the highresolution 3D image data. These simulations assist in the characterization of the effective mechanical behavior.
Due to the sheer size of today’s CT images, the resulting computational homogenization algorithms face severe challenges related to the high computational demands. For instance, to compute effective linear elastic properties on a 4096^{3}
CT image, the number of nodal displacement degrees of freedom amounts to approximately 206 Billion. The solution of problems of this size using conventional simulation methods such as the finite element (FE) method requires huge compute clusters
Arbenz2008; Arbenz2014. These difficulties are commonly overcome by working with conventional FEM on a variety of smaller subsamples of moderate size Kanit:2003. The resulting effective properties of the individual subsamples are averaged in postprocessing, cf., for instance, Andrae2013a; Andrae2013b. This method, however, does not exploit the available information gathered from the specimen.FFTbased solvers
In contrast to the conventional FEM approach, the numerical homogenization method of MoulinecSuquet MoulinecSuquet1994; MoulinecSuquet1998
operates on the voxels of a CT image directly: the set of unknowns is formed by the strains, i.e. one tensor for each voxel in the CT image. The solution algorithm works in place (i.e. matrixfree) such that, in practice, the size of the CT images to be treated is only restricted by the size of the memory and the affordable compute time. FFTbased schemes can handle arbitrary heterogeneity, phase contrast
(e.g., MichelMoulinecSuquet2001)and degree of compressibility of the materials. Generally, the number of iterations is independent of the grid size, depending only on the material’s contrast, i.e., the maximum of the quotient of the largest and the smallest eigenvalue of the elastic tensor field, and on the geometric complexity. The solution of the linear algebraic system relies upon a LippmannSchwinger fixed point formulation, enabling, by use of fast Fourier transform (FFT), a fast matrixfree implementation. By changing the discretization to finite elements
Willot2015 and finite differences Schneider2016 also infinite material contrast problems became solvable. Recent displacementbased implementations Leuschner2018 also allowed the use of higherorder integration schemes without increasing the memory demands and with improved computational efficiency over strainbased algorithms also building on finite element technology Schneider2017; Leuschner2018. A relation with Galerkin based methods was outlined by Vondrejc:2014 where the use of the efficient Conjugate Gradient (CG) method was mentioned. Interestingly, the CG method (as well as minres and other krylov solvers) is completely natural in FANS and FFT Hex Leuschner2018; Schneider2017, including straightforward implementation.Composite voxel technique
Despite the many advances regarding theory and algorithmic realization of FFTbased methods, the resolution of stateoftheart computed tomography poses challenges for the numerics. A single, doubleprecision scalar field on a 4096^{3} voxel image as delivered by modern CT scanners already occupies 512 GB of memory. Thus, even performing linear elastic computations on such images requires either exceptionally equipped workstations or big compute clusters. This applies evermore so to inelastic computations, where additional history variables need to be stored, increasing the memory demand significantly.
To enable computations on conventional desktop computers or workstations that can still take into account relevant microstructural details of the fully resolved image, the composite voxel technique was developed for linear elastic, hyperelastic and inelastic problems LionelGelebart2015; Kabel2015; Kabel2016; Kabel2017. A coarsegraining procedure serves as the initial idea, i.e., a number of smaller, typically , voxels are merged into bigger voxels. Each of these socalled composite voxels
gets assigned an appropriately chosen material law which is based on a biphasic laminate. This laminate takes into account the exact phase volume fractions and the interface normal vector
. Thereby, it can reflect the most relevant microstructural features, avoid staircase phenomena intrinsic to regular voxel discretizations and allow for a boost in accuracy in computational homogenization.In many scenarios the edge lengths of the voxels are different. This can be due to the employed imaging technique, e.g., in serial sectioning by a focused ion beam and using scanning electron microscopy (FIBSEM, Uchic:2006). Another reason for the wish of using anisotropic voxel shapes is the presence of a preferred direction, e.g., in composites with preferred orientation such as long fiber reinforced materials studied by Fliegener:2014. The consideration of anisotropic composite voxels – called boxels in the present study – leads to low quality normal orientations when using the original approach suggested by Kabel:2015. Likewise, low volume fractions can deteriorate the quality of the normal orientation.
Another issue arises depending on the material contrast and the local volume fraction in the composite voxels: If the soft phase of the twophase laminate has a small volume fraction, the deformation within this phase is severely exagerated. This typically leads to convergence issues of the NewtonRaphson method which evaluates the effective response of the composite voxel and is normally circumvented by limiting the volume fractions of either of constituents to be in the range of 595% and by resorting to simple Voigt averaging otherwise, ruling out the aforementioned advantages of the compositve voxels.
Outline
In the current study we target composite voxels for finite strain homogenization problems. In Section 2 the homogenization problem is recalled. The foundations of FFTbased schemes are summarized in Section 3 including the LippmannSchwinger equation and an extension of the staggered grid approach of Schneider2016 towards improved local fields. The composite voxel technique is considered in detail in Section 4 including algorithmic improvements leading to increased robustness for finite strain problems.
Notation
The spatial average of a quantity over a domain with measure is defined as
(1) 
In the sequel bold face lowercase letters denote vectors (exceptions: normal vector and traction ), bold face upper case letters denote 2tensors and blackboard bold uppercase letters (e.g., ) denote 4tensors. The inner product contracts all indices of two tensors while the usual linear mapping contracts the last indices of the first tensor with the following tensor. The double contraction operator contracts the trailing two indices of with the leading two indices of . The outer tensor product is denoted by and its symmetric version is given by .
General vectors and matrices are denoted by single and double underlines, respectively. Note that in the sequel we focus on an orthonormal basis for the tensor fields. A vector representation for general () and symmetric 2tensors () is given by
(2) 
Note the different ordering and normalization factors of the Mandeltype notation versus the commonly used Voigt notation which is avoided due to the ambiguity of stress vs. strain like quantities in the vector representation. The notation (2) preserves the inner product for general 2tensors (e.g., ) and symmetric 2tensors (e.g., ), i.e.
(3) 
Likewise, 4tensors are represented by and (in the case of left and right subsymmetry) matrices, respectively.
In view of spatial derivatives, the gradient () and the divergence () with respect to the reference configuration are used. In the context of an internal surface the jump operator is given as
(4) 
where is the normal of pointing from the second phase (referred to as inclusion phase) into the first phase (referred to as the matrix material).
2 Homogenization Problem
We focus our attention to finite strain mechanics of solids to describe the framework. Additinally, we explain the simplifications for infinitesimal strain mechanics. Consider a material point of a body in the reference configuration at time . The corresponding current position at time is denoted by in where is the deformed domain. The motion of the body is given by
(5) 
where is the spatial dimension. In the sequel, the dependence on time is implicitly assumed but not reflected in the arguments for brevity. The displacement field is given by . The displacement gradient and the deformation gradient are defined as
(6)  
(7) 
The volume change at a material point is given by the material Jacobian of the deformation map,
(8) 
with the differential volume in the current configuration and is its counterpart in the reference configuration. The right CauchyGreen tensor is given by . At a particular material point with an infinitesimal area with a unit normal vector , we define the resulting traction vector in terms of the first PiolaKirchoff (PK1) stress tensor ,
(9) 
where the first PiolaKirchoff (PK1) stress can be expressed in terms of the deformation gradient and the second PiolaKirchoff stress (PK2) via .
In homogenization, the constitutive response on the structural (or macroscopic) scale of a body depends strongly on the underlying microstructure. In the following we consider a statistically homogeneous periodic microstructure that is described by a periodic representative volume element (RVE)
(10) 
on the microscale. Separation of the characteristic microscopic and macroscopic lengths and , respectively, is assumed. Formally, this is equivalent to requiring
(11) 
In the absence of scale separation, advanced modeling techniques must be taken into account that are beyond the scope of the current investigation, e.g., by using filters Yvonnet2014 or higherorder continuum theories (e.g., Jaenicke2009).
The objective of homogenization is the identification of the effective, macroscopic constitutive response of microheterogeneous materials in order to render twoscale simulations feasible. Equilibrium on the microstructural RVE induces the virtual work principle. The latter implies the crucial HillMandel macrohomogeneity condition:
(12) 
In the infinitesimal strain framework this is equivalent to
(13) 
where is the Cauchy stress tensor and is the infinitesimal strain tensor.
Periodic fluctuation boundary conditions which satisfy the HillMandel condition Suquet1985 are used. We decompose the boundary of the RVE into two parts and , such that , where the material points corresponding to matching points on and are denoted by and , respectively. For an imposed macroscopic deformation over the RVE the boundary conditions induce
(14)  
(15) 
for point pairs . The displacement can be decomposed into a rigid translation of the RVE, a homogeneous field and a displacement fluctuation . Condition (14) induces the periodicity of the fluctuations, i.e.
(16) 
The related function space for is referred to as which, in the case of finite element simulations, is a subset of the usual Sobolev space of weakly differentiable functions . The rigid translation does not influence the homogenization and can be set to . In a Lagrangian setting with periodic boundary conditions and (static) equilibrium condition without volume forces, the homogenization problem (HOM) on for prescribed deformation reads:
(17)  
(18)  
(19)  
(20) 
The soughtafter output of (HOM) is the homogenized macroscopic stress
(21) 
3 FFTbased homogenization
3.1 Foundations for FFTbased schemes
The solution of (HOM) (17)(20) can either be obtained (in seldom cases) using analytical solution or by using discrete numerical techniques. The most significant methods for computational homogenization in solid mechanics are certainly the finite element method (FEM) (e.g., FE^{2}, Feyel1999) and FFTbased homogenization.
The latter was proposed by Moulinec and Suqet MoulinecSuquet1994 in the middle of the 1990s for problems of linear elasticity. It is based on the LippmannSchwinger equation in elasticity Kroener1977; ZellerDederichs1973 and avoids both, time consuming meshing needed by conforming finite element discretizations as well as the assembly of the related linear system. Therefore, the memory needed for solving the problem is significantly reduced compared with other methods. By virtue of the seminal Fast Fourier Transform (FFT, Cooley1965) algorithm, the compute time scales with where is the number of unknowns, i.e. it is just slightly superlinear.
In view of large deformations Lahellec, Moulinec, and Suquet Lahellec2001 proposed to solve the nonlinear LippmannSchwinger equation for finite strains by the NewtonRaphson method and the linear MoulinecSuquet fixed point solver. In contrast, Eisenlohr et al. Eisenlohr2013 suggested using the MoulinecSuquet fixed point iteration on the nonlinear LippmannSchwinger equation for finite strains directly. Kabel et al. Kabel2014 carried over the idea of Vinogradov and Milton VinogradovMilton2008 and of Gélébart and MondonCancel GelebartMondonCancel2013 to combine the NewtonRaphson procedure with fast linear solvers to the geometrically nonlinear case.
The present work will also make use of the nonlinear conjugate gradient (CG) method introduced by Schneider Schneider2020 for small strains. Further, the FFTaccelerated solution of regular trilinear hexahedral elements (FFT Hex cf., Schneider2017) and in the closely related FourierAccelerated Nodal Solver Leuschner2018 will be used in the sequel, for which the nonlinear CG method has a perfectly natural interpretation as a preconditioner that is more subtle for other FFT schemes. An overview on alternative solution methods can be found in Schneider2021.
3.2 The LippmannSchwinger equation in hyperelasticity
The core of FFTbased homogenization schemes is the LippmannSchwinger equation in hyperelasticity which are briefly summarized in the following.
Let a constant reference stiffness satisfying
(22) 
and a positive constant independent of . Introducing the stress polarization
(23) 
the equilibrium equation (18) transforms into the relation
(24) 
In the sequel is the Sobolev space of periodic functions with zero mean on that are weakly differentiable and is its dual. Then we denote by the solution operator of the linear reference problem on , which associates to a righthand side the solution of the variational equation
(25) 
where . Then
(26)  
(27) 
or, equivalently,
(28) 
with the bounded linear operator defined by
(29) 
which is usually called Green’s operator. The LippmannSchwinger equation (28) was first introduced in Kroener1971; ZellerDederichs1973 and is equivalent to the variational form Vondrejc2013
(30) 
of the equilibrium equation.
3.3 A finite difference discretization on a staggered grid
In the sequel, we describe the consistent formulation of the staggered grid discretization Schneider2016 for the geometrically nonlinear case. Similar to the small strain case, the diagonal and offdiagonal components of the deformation gradient are located at different positions, compare Fig. 1. It is shown that ignoring this leads to unsatisfactory results which can be significantly improved by using a doublyfine material grid (DFMG). In contrast to the linear elastic small strain regime, however, the DFMG needs an increased number of evaluations of the nonlinear material law.
Fix positive integers and consider a regular periodic grid consisting of cells, each a translate of with (). Let and denote the location of the staggered grid coordinates, i.e. for . For the cell the coordinates of the displacement vector are located at the cell face centers, i.e. is located at the coordinate , lives at and is associated to . The situation is displayed in Figure 1 (a), where we have restricted ourselves to the D case for clarity.
The diagonal components () of the deformation gradient are positioned on the cell centers , whereas the offdiagonal strains , and , are located on the corresponding edge midpoints , and . For visualization, we again refer to Figure 1 (a).
Displacements and gradients are connected by central difference formulae, where periodicity is understood implicitly. More precisely, introduce forward and backward difference operators on a scalar discrete field by the formulae
(31) 
for . Then we introduce the gradient operator
(32) 
giving rise to the deformation gradient associated to a periodic displacement field
Similarly, there is a divergence operator, turning into
(33) 
The stress variable is located on the same position as the corresponding .
The constitutive law is defined piecewise on the regular voxel grid. For a typical material cell, shaded in gray in Figure 1 (b), it is not clear how to apply the material law for offdiagonal strains, as the function may be defined differently along the corresponding edge.
We circumvent these problems by utilizing a doublyfine grid, i.e. a grid with half the spacing of the original grid. Figure 1 (b) illustrates this concept – a typical doubly fine cell is shaded in orange. We interpret the deformation gradients and stresses as living on this doublyfine grid. For every deformation gradient component and every doublyfine cell there is precisely one value, as specified in Figure 1 (a), located on the boundary of the cell. We associate this value to the doublyfine cell. Thus, a particular value is distributed to the (in D) or (in D) adjacent doublyfine cells, compare Figure 2 for an illustration. The stress components are distributed similarly to the deformation gradients, i.e. the staggering of Figure 2, is also present for the stresses.
With this assignment, to any doublyfine grid cell all (in D) or (in D) deformation gradient components are associated. The discretization just outlined directly carries over to three space dimensions easily. The variable placement in this case is shown in Figure 3.
We suppose that the constitutive material law is given on the original grid, i.e. each cell is associated with one material. We will elaborate on the implementation of the material law. By averaging over the combination of adjacent doublyfine voxels, the material law can be written as
(34) 
with
(35)  
(36) 
and
for
as well as
where
4 Overcoming regular grids: Composites Voxels
One of the main advantages but  at the same time  a major disadvantage of FFTbased techniques is the constraint to regular Cartesian grids: They are unable to capture the material interface accurately compared to interfaceconforming discretizations. This is primarily due to the binarized nature of the microstructure which leads to the so called
staircase approximation of the interface. In order to capture the microscale effects sufficiently, a high grid resolution is needed which–in turn–calls for high computational cost. In order to eliminate or at least limit the staircase phenomenon without the need of a (distinct) grid refinement, socalled composite voxels were previously suggested in Kabel2015; Merkert2015; Kabel2016; Kabel2017. They enhance the existing binary discretization with special interface voxels with effective material properties that depend on the phase volume fractions and the surface normal through a laminate problem.Consider a interface voxel comprised of two phases. The two phases are denoted by and , respectively, where and stand for inclusion and matrix phase, respectively. The fields corresponding to the two phases are represented as for brevity and either phase is assumed to be equipped with a hyperelastic strain energy density . Within material interface is approximated with a plane leading to a rank1 laminate defined by the interface normal vector in material configuration. The individual phase volume fractions are given by and where .
In previous works, different rules of mixture have been suggested, namely the Voigt (•^{V}) and Reuss (•^{R})^{*}^{*}*In the nonlinear kinematic regime these correspond to the Taylor and Sachs approximation, respectively.estimates corresponding to upper and lower bounds of the mechanical response. Hill’s estimate (•^{H}) corresponds to the arithemetic mean of the stiffness of the Voigt and Reuss estimate:
(37)  
(38)  
(39) 
Obviously the Voigt/Taylor upper and Reuss/Sachs lower bounds are overly stiff and soft, respectively, due to simplistic physical considerations, i.e., the assertion of constant deformation and stress, respectively. Hill’s estimate is a practical compromise with limited physical backing.
Inspired by laminate theory, Milton2002 suggests the effective elasticity tensor of the laminate to be defined implicitly by
(40) 
where the fourth order tensor depends on the normal according to
(41) 
for . Note that in laminate theory the states in the phase are piecewise constant. Hence, the averaging translates into
(42) 
The laminate mixing rule takes the interface orientation into consideration and introduces a parameter devoid of immediate physical interpretation that must be chosen carefully. More precisely, has to be larger than the largest eigenvalue of the individual stiffness tensors . Solving the laminate mixing rule (4) for linear elastic materials involves inversions of matrices^{†}^{†}†Symmetry is exploited.. Note that this effort is needed for each composite voxel individually, i.e., at every interface related voxel or cubature point within an FFTbased simulation, which can lead to nonnegligible overhead. Hence, a more numerically efficient and physics motivated approach devoid of additional parameters is sketched in the following that can also handle geometric and material nonlinearities, if needed.
The composite voxel concept has also been implemented for steady state heat conduction problems, too. In the interest of compactness we provide all that is needed for their implementation in Appendix 5.
4.1 Hadamard jump conditions for finite strain kinematics
Following established laminate theory, the following assumptions hold:

The fields are constant in each subdomain

The deformation gradient in and are related by a rank1 jump along the normal orientation, i.e. for ,
(43) 
The traction vector is continuous across the interface , i.e.,
(44) (45)
The kinematic compatibility on the interface is, thus, characterized by having a continuous deformation gradient in the tangential direction to the interface. This ensures that neither material penetration nor cracking occurs. Condition (45) ensures that the interface is in static equilibrium.
In order to enforce (43), we chose the following parameterization of the deformation gradient tensors of the two material phases
(46) 
where is related to in (43) by a scaling constant and is the prescribed deformation gradient on the composite voxel . The parameterization (46) preserves the volume average of the deformation gradient over the composite voxel according to
(47) 
In the finite strain context not just the average of must be preserved, but also the total volume must remain constant. This is stated by the following Lemma.
Lemma 1.
The rank1 kinematic perturbation following the Hadamard jump condition on the material interface is volume preserving by construction.
Proof.
Suppose is an invertible square matrix and , are column vectors. Then the matrix determinant lemma states that
(48) 
The deformation gradient is regular by definition since material inversion is not allowed. Thus, by setting and the material Jacobian in the two phases can be computed using the matrix determinant lemma
(49) 
Averaging over the composite voxel , we obtain
(50) 
∎
4.2 Motivation via linearized kinematics
In the context of small deformations, the infinitesimal strain tensor
(53) 
and the Cauchy stress are used. In this setting and in the case of linear elasticity, (51) can be solved analytically to obtain closed form solutions. Balancing the tractions on yields
(54)  
(55) 
The operator is the symmetric part of the dyadic product. Switching to a matrixvector notation for symmetric tensors and expressing in Mandel notation^{‡}^{‡}‡(with respect to the following orthonormal basis of the space of symmetric tensors: , see Notation section after the Introduction), the perturbation strain due to the rank1 jump is given by (here without loss of generality)
(56) 
Substitution of (55) into (54) yields equations for unknowns namely the components of the vector defining the strain perturbation within the laminate. Solving for we get
(57)  
(58) 
Making use of the explicit solution (57) the effective stress tensor
(59) 
is gained. By expanding
(60) 
we obtain the closed form analytical solution of a laminate material constitutive tangent moduli for a linear elastic (but generally anisotropic) biphasic laminate:
(61)  
(62) 
The closed form solution in (62) requires just a single matrix inversion per composite voxel for a linear material model. It is devoid of additional parameters and, technically, it can be precomputed offline for all composite voxels. Note the straightforward symmetric structure of which consists of:

the Voigt stiffness,

the difference of the material stiffness tensors,

and a matrix that can be interpreted as a compliance.
This formulation can also be easily extended to a purely thermal problem with linear thermal conductivity tensors, see Appendix 5 for details.
4.2.1 Generalization to nonlinear kinematics
In the case of large strain kinematics, the system (51) gets nonlinear and, hence, it must be solved iteratively using, e.g., a NewtonRaphson scheme with the Hessian as
(63) 
The Hessian in index notation can be written in terms of the 2tensors and as