FFT-based Homogenization at Finite Strains using Composite Boxels (ComBo)

by   Sanath Keshav, et al.

Computational homogenization is the gold standard for concurrent multi-scale simulations (e.g., FE2) in scale-bridging applications. Experimental and synthetic material microstructures are often represented by 3D image data. The computational complexity of simulations operating on such three-dimensional high-resolution voxel data comprising billions of unknowns induces the need for algorithmically and numerically efficient solvers. The inability of voxelized 3D geometries to capture smooth material interfaces accurately, along with the necessity for complexity reduction, motivates a special local coarse-graining technique called composite voxels [Kabel,M. et al. (2015)]. Composite voxels condense multiple fine-scale voxels into a single voxel obeying a theory-inspired constitutive model by employing laminate theory. Composite voxels enhance local field quality at a modest computational cost. Our contribution comprises the generalization towards composite boxels (ComBo) that are nonequiaxed, a feature that can pay off for materials with a preferred direction. A novel image-based normal detection algorithm is devised which improves the accuracy by around 30% against the orientation cf. [Kabel,M. et al. (2015) ]. Further, the use of ComBo for finite strain simulations is studied in detail. An efficient implementation is proposed, and an essential back-projection algorithm preventing physically inadmissible states is developed, which improves robustness. Various examples show the efficiency of ComBo and the proposed algorithmic enhancements for nonlinear mechanical problems. The general usability is emphasized by examining and comparing the performance of myriad Fast Fourier Transform (FFT) based solvers including a detailed description of the new Doubly-Fine Material Grid (DFMG). All of the employed schemes benefit from the ComBo discretization.


page 14

page 15

page 18

page 20

page 22

page 23

page 25

page 26


An FE-DMN method for the multiscale analysis of thermomechanical composites

We extend the FE-DMN method to fully coupled thermomechanical two-scale ...

Thermoplastic Fiber Reinforced Composite Material Characterization and Precise Finite Element Analysis for 4D Printing

Four-dimensional (4D) printing, a new technology emerged from additive m...

Robust concurrent topology optimization of structure and its composite material considering uncertainty with imprecise probability

This paper studied a robust concurrent topology optimization (RCTO) appr...

QFT-based Homogenization

Efficient numerical characterization is a key problem in composite mater...

A geometrically adapted reduced set of frequencies for a FFT-based microstructure simulation

We present a modified model order reduction (MOR) technique for the FFT-...

Multi-physics integration platform MuPIF: Application for composite material design

This paper presents the design of the MuPIF distributed, multi-physics s...

1 Introduction

Homogenization in engineering

In the last decade, the quality of micro x-ray 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 40963 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 micro-cracks. In order to understand the effect of such microscopic features, one has to solve PDEs on the high-resolution 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 40963

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 post-processing, cf., for instance, Andrae2013a; Andrae2013b. This method, however, does not exploit the available information gathered from the specimen.

FFT-based solvers

In contrast to the conventional FEM approach, the numerical homogenization method of Moulinec-Suquet 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. matrix-free) 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. FFT-based 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 Lippmann-Schwinger fixed point formulation, enabling, by use of fast Fourier transform (FFT), a fast matrix-free implementation. By changing the discretization to finite elements

Willot2015 and finite differences Schneider2016 also infinite material contrast problems became solvable. Recent displacement-based implementations Leuschner2018 also allowed the use of higher-order integration schemes without increasing the memory demands and with improved computational efficiency over strain-based 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 straight-forward implementation.

Composite voxel technique

Despite the many advances regarding theory and algorithmic realization of FFT-based methods, the resolution of state-of-the-art computed tomography poses challenges for the numerics. A single, double-precision scalar field on a 40963 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 coarse-graining procedure serves as the initial idea, i.e., a number of smaller, typically , voxels are merged into bigger voxels. Each of these so-called composite voxels

gets assigned an appropriately chosen material law which is based on a bi-phasic 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 (FIB-SEM, 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 two-phase laminate has a small volume fraction, the deformation within this phase is severely exagerated. This typically leads to convergence issues of the Newton-Raphson 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.


In the current study we target composite voxels for finite strain homogenization problems. In Section 2 the homogenization problem is recalled. The foundations of FFT-based schemes are summarized in Section 3 including the Lippmann-Schwinger 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.


The spatial average of a quantity over a domain with measure is defined as


In the sequel bold face lowercase letters denote vectors (exceptions: normal vector and traction ), bold face upper case letters denote 2-tensors and blackboard bold uppercase letters (e.g., ) denote 4-tensors. 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 2-tensors () is given by


Note the different ordering and normalization factors of the Mandel-type 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 2-tensors (e.g., ) and symmetric 2-tensors (e.g., ), i.e.


Likewise, 4-tensors are represented by and (in the case of left and right sub-symmetry) 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


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


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


The volume change at a material point is given by the material Jacobian of the deformation map,


with the differential volume in the current configuration and is its counterpart in the reference configuration. The right Cauchy-Green 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 Piola-Kirchoff (PK1) stress tensor ,


where the first Piola-Kirchoff (PK1) stress can be expressed in terms of the deformation gradient and the second Piola-Kirchoff 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)


on the microscale. Separation of the characteristic microscopic and macroscopic lengths and , respectively, is assumed. Formally, this is equivalent to requiring


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 higher-order continuum theories (e.g., Jaenicke2009).

The objective of homogenization is the identification of the effective, macroscopic constitutive response of micro-heterogeneous materials in order to render twoscale simulations feasible. Equilibrium on the microstructural RVE   induces the virtual work principle. The latter implies the crucial Hill-Mandel macro-homogeneity condition:


In the infinitesimal strain framework this is equivalent to


where is the Cauchy stress tensor and is the infinitesimal strain tensor.

Periodic fluctuation boundary conditions which satisfy the Hill-Mandel 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


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.


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:


The sought-after output of (HOM) is the homogenized macroscopic stress


3 FFT-based homogenization

3.1 Foundations for FFT-based 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., FE2, Feyel1999) and FFT-based 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 Lippmann-Schwinger 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 Lippmann-Schwinger equation for finite strains by the Newton-Raphson method and the linear Moulinec-Suquet fixed point solver. In contrast, Eisenlohr et al. Eisenlohr2013 suggested using the Moulinec-Suquet fixed point iteration on the nonlinear Lippmann-Schwinger equation for finite strains directly. Kabel et al. Kabel2014 carried over the idea of Vinogradov and Milton VinogradovMilton2008 and of Gélébart and Mondon-Cancel GelebartMondonCancel2013 to combine the Newton-Raphson 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 FFT-accelerated solution of regular tri-linear hexahedral elements (FFT- Hex cf., Schneider2017) and in the closely related Fourier-Accelerated 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 Lippmann-Schwinger equation in hyperelasticity

The core of FFT-based homogenization schemes is the Lippmann-Schwinger equation in hyperelasticity which are briefly summarized in the following.

Let a constant reference stiffness satisfying


and a positive constant independent of . Introducing the stress polarization


the equilibrium equation (18) transforms into the relation


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 right-hand side the solution of the variational equation


where . Then


or, equivalently,


with the bounded linear operator defined by


which is usually called Green’s operator. The Lippmann-Schwinger equation (28) was first introduced in Kroener1971; ZellerDederichs1973 and is equivalent to the variational form Vondrejc2013


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 off-diagonal 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 doubly-fine 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.

Figure 1: Placement of the strain and displacement variables in D: (a) Location of strain and displacement variables. (b) The variables’ grid (gray box) and the doubly fine material grid (orange box).

The diagonal components () of the deformation gradient are positioned on the cell centers , whereas the off-diagonal 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


for . Then we introduce the gradient operator


giving rise to the deformation gradient associated to a periodic displacement field

Similarly, there is a divergence operator, turning into


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 off-diagonal strains, as the function may be defined differently along the corresponding edge.

Figure 2: Placement of deformation gradient variables within the doubly-fine material grid: (a) Location of (shaded in orange). (b) Location of (shaded in orange). The variables’ grid is depicted by the dashed black lines.

We circumvent these problems by utilizing a doubly-fine 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 doubly-fine grid. For every deformation gradient component and every doubly-fine 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 doubly-fine cell. Thus, a particular value is distributed to the (in D) or (in D) adjacent doubly-fine 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 doubly-fine 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.

Figure 3: Placement of the deformation gradients and displacement variables in D.

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 doubly-fine voxels, the material law can be written as






as well as


4 Overcoming regular grids: Composites Voxels

One of the main advantages but - at the same time - a major disadvantage of FFT-based techniques is the constraint to regular Cartesian grids: They are unable to capture the material interface accurately compared to interface-conforming 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, so-called 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 rank-1 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:


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


where the fourth order tensor depends on the normal  according to


for . Note that in laminate theory the states in the phase are piecewise constant. Hence, the averaging translates into


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 matricesSymmetry 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 FFT-based simulation, which can lead to non-negligible 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 sub-domain

  • The deformation gradient in and are related by a rank-1 jump along the normal orientation, i.e. for ,

  • The traction vector is continuous across the interface , i.e.,


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


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


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 rank-1 kinematic perturbation following the Hadamard jump condition on the material interface is volume preserving by construction.


Suppose is an invertible square matrix and , are column vectors. Then the matrix determinant lemma states that


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


Averaging over the composite voxel , we obtain


To obtain the gradient jump vector the equilibrium of the tractions on the interface cf. (45) must be granted:


Examining the relative volume in , additional constraints on emerge:


The two constraints due to (52) will gain in importance at a later stage.

4.2 Motivation via linearized kinematics

In the context of small deformations, the infinitesimal strain tensor


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


The operator is the symmetric part of the dyadic product. Switching to a matrix-vector 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 rank-1 jump is given by (here without loss of generality)


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


Making use of the explicit solution (57) the effective stress tensor


is gained. By expanding


we obtain the closed form analytical solution of a laminate material constitutive tangent moduli for a linear elastic (but generally anisotropic) biphasic laminate:


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 straight-forward 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 Newton-Raphson scheme with the Hessian  as


The Hessian in index notation can be written in terms of the 2-tensors and as