Computational Mechanics (2023) 71:191–212 https://doi.org/10.1007/s00466-022-02232-4 ORIG INAL PAPER FFT-based homogenization at finite strains using composite boxels (ComBo) Sanath Keshav1 · Felix Fritzen1 ·Matthias Kabel2 Received: 28 April 2022 / Accepted: 15 September 2022 / Published online: 10 October 2022 © The Author(s) 2022 Abstract Computational homogenization is the gold standard for concurrent multi-scale simulations (e.g., FE2) in scale-bridging applications. Often the simulations are based on experimental and synthetic material microstructures represented by high- resolution 3D image data. The computational complexity of simulations operating on such voxel data is distinct. The inability of voxelized 3Dgeometries to capture smoothmaterial interfaces accurately, alongwith the necessity for complexity reduction, has motivated a special local coarse-graining technique called composite voxels (Kabel et al. Comput Methods Appl Mech Eng 294: 168–188, 2015). They condense multiple fine-scale voxels into a single voxel, whose constitutive model is derived from the laminate theory. Our contribution generalizes composite voxels towards composite boxels (ComBo) that are non- equiaxed, a feature that can pay off for materials with a preferred direction such as pseudo-uni-directional fiber composites. A novel image-based normal detection algorithm is devised which (i) allows for boxels in the firsts place and (ii) reduces the error in the phase-averaged stresses by around 30% against the orientation cf. Kabel et al. (Comput Methods Appl Mech Eng 294: 168–188, 2015) even for equiaxed voxels. Further, the use of ComBo for finite strain simulations is studied in detail. An efficient and robust implementation is proposed, featuring an essential selective back-projection algorithm preventing physically inadmissible states. Various examples show the efficiency of ComBo against the original proposal by Kabel et al. (ComputMethods ApplMech Eng 294: 168–188, 2015) and the proposed algorithmic enhancements for nonlinear mechanical problems. The general usability is emphasized by examining various Fast Fourier Transform (FFT) based solvers, including a detailed description of the Doubly-Fine Material Grid (DFMG) for finite strains. All of the studied schemes benefit from the ComBo discretization. Keywords Homogenization · Fast fourier transform · Composite voxel · Composite boxel B Felix Fritzen fritzen@simtech.uni-stuttgart.de Sanath Keshav keshav@mib.uni-stuttgart.de Matthias Kabel matthias.kabel@itwm.fraunhofer.de 1 SC Simtech, Data Analytics in Engineering, University of Stuttgart, Universitätsstr.32, 70569 Stuttgart, Germany 2 Department Flow and Material Simulation, Fraunhofer-Institut für Techno- und Wirtschaftsmathematik ITWM, 67663 Kaiserslautern, Germany 1 Introduction 1.1 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 compos- ite 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- 123 http://crossmark.crossref.org/dialog/?doi=10.1007/s00466-022-02232-4&domain=pdf http://orcid.org/0000-0002-9131-4560 http://orcid.org/0000-0003-4926-0068 http://orcid.org/0000-0001-8869-2387 192 Computational Mechanics (2023) 71:191–212 resolution 3D image data. These simulations assist in the characterization of effective mechanical behavior. Due to the sheer size of today’s CT images, the result- ing 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 solu- tion to problems of this size using conventional simulation methods such as the finite element (FE) method requires huge compute clusters [1,2]. These difficulties are commonly overcome by working with conventional FEM on a variety of smaller subsamples of moderate size [3]. The resulting effec- tive properties of the individual subsamples are averaged in post-processing, cf., for instance, [4,5]. This method, how- ever, does not exploit the available information gathered from the specimen. 1.2 FFT-based solvers In contrast to the conventional FEM approach, the numerical homogenization method of Moulinec–Suquet [6,7] 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., [8]], 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 trans- form (FFT), a fast matrix-free implementation. By changing the discretization to finite elements [9] and finite differ- ences [10] also infinite material contrast problems became solvable. Recent displacement-based implementations [11] 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 onfinite element technology [11,12].A relationwith Galerkin-based methods was outlined by [13] where the use of the efficient Conjugate Gradient (CG) method was men- tioned. Interestingly, the CG method (as well as minres and other Krylov solvers) is completely natural in FANS, and FFT-Q1 Hex [11,12], including straightforward implemen- tation. 1.3 Composite voxel technique Despite themany advances in theory and algorithmic realiza- tion of FFT-based methods, the resolution of state-of-the-art computed tomography poses challenges for numerics. A sin- gle, double-precision scalar field on a 40963 voxel image as delivered bymodern µCT scanners already occupies 512 GB of memory. Thus, even performing linear elastic computa- tions on such images requires either exceptionally equipped workstations or big compute clusters. This applies evermore so to inelastic computations, where additional history vari- ables need to be stored, increasing the memory demand significantly. To enable computations on conventional desktop comput- ers or workstations that can still take into account relevant microstructural details of the fully resolved image, the com- posite voxel technique was developed for linear elastic, hyperelastic and inelastic problems [14–17]. A coarse- graining procedure serves as the initial idea, i.e., a number of smaller, typically 43 = 64, voxels are merged into big- ger voxels. Each of these so-called composite voxels gets assigned an appropriately chosenmaterial lawwhich is based on a bi-phasic laminate. This laminate takes into account the exact phase volume fractions and the interface normal vec- tor N . Thereby, it can reflect themost relevantmicrostructural features, avoid staircase phenomena intrinsic to regular voxel discretizations and allow for a boost in accuracy in compu- tational homogenization. In many scenarios, the edge lengths of the voxels are dif- ferent. 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, [18]). Another reason for the wish to use anisotropic voxel shapes is the presence of a preferred direction, e.g., in composites with a preferred orientation such as long fiber reinforced materials studied by [19]. 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 [15]. Likewise, low volume fractions can dete- riorate 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 vol- ume fraction, the deformation within this phase is severely exaggerated. This typically leads to convergence issues of the Newton–Raphson method, which evaluates the effective response of the composite voxel and is normally circum- vented by limiting the volume fractions of either constituent to be in the range of 5· · · 95%and by resorting to simpleVoigt averaging otherwise, ruling out the aforementioned advan- tages of the composite voxels. 123 Computational Mechanics (2023) 71:191–212 193 1.4 Outline In the current study, we target composite voxels for finite strain homogenization problems. In Sect. 2 the homoge- nization problem is recalled. The foundations of FFT-based schemes are summarized in Sect. 3 including the Lippmann– Schwinger equation and an extension of the staggered grid approach of [10] towards improved local fields. The compos- ite voxel technique is considered in detail in Sect. 4 including algorithmic improvements leading to increased robustness for finite strain problems. 1.5 Notation The spatial average of a quantity over a domainAwith mea- sure A = |A| is defined as 〈·〉A = 1 A ∫ A (·) dA . (1) In the sequel boldface lowercase letters denote vectors (exceptions: material coordinate X , normal vector N , dis- placement U and traction T ), boldface upper case letters denote 2-tensors and blackboard bold uppercase letters (e.g., C) denote 4-tensors. The inner product contracts all indices of two k-tensors while the usual linear mapping contracts the last k indices of the first tensor with the following k-tensor. The double contraction operator A : B contracts the trailing two indices ofAwith the leading two indices of B. The outer tensor product is denoted by ⊗, and its symmetric version is given by ⊗s. 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 (F ↔ F) and symmetric 2-tensors (S ↔ S) is given by F = ⎡ ⎢⎢⎢⎢⎢⎢⎢⎣ F11 F12 F13 F21 ... F33 ⎤ ⎥⎥⎥⎥⎥⎥⎥⎦ ∈ R 9, S = ⎡ ⎢⎢⎢⎢⎢⎢⎣ S11 S22 S33√ 2S12√ 2S13√ 2S23 ⎤ ⎥⎥⎥⎥⎥⎥⎦ ∈ R 6. (2) Note the different ordering and normalization factors of the Mandel-type notation versus the commonly used Voigt nota- tion, which is avoided due to the ambiguity of stress versus strain-like quantities in the vector representation. The nota- tion (2) preserves the inner product for general 2-tensors (e.g., F, P) and symmetric 2-tensors (e.g., S, E), i.e., F · P = FTP, S · E = STE . (3) Likewise, 4-tensors are represented by 9×9 and 6×6 (in the case of left and right sub-symmetry) matrices, respectively. In view of spatial derivatives, the gradient (Grad (·)) and the divergence (Div (·)) with respect to the reference config- uration are used. In the context of an internal surfaceS , the jump operator is given as �u� = lim ε→0+ u(x + εN) − u(x − εN) (x ∈ S ), (4) where N is the normal of S pointing from the second phase (referred to as the inclusion phase) into the first phase (referred to as the matrix material). 2 Homogenization problem Wefocus our attention on thefinite strainmechanics of hyper- elastic solids. Consider a material point X of a body B0 in the reference configuration at time t = 0. The corresponding current position at time t ∈ [0, T ] is denoted by x in the deformed domain B(t). The motion of the body is given by ϕ(X, t) : B0 × R → R d , x = ϕ(X, t), (5) where d ∈ {2, 3} is the spatial dimension. In the sequel, the dependence on time t is implicitly assumed but not reflected in the arguments for brevity. The displacement field is given by U(X) = x − X . The deformation gradient F is defined as F = Grad (ϕ(X)) = I + H . (6) The volume change at amaterial point is given by thematerial Jacobian of the deformation map, dv dV = J = det F > 0, (7) with dv the differential volume in the current configuration and dV denoting its counterpart in the reference configura- tion. The right Cauchy-Green tensor is given by C = FTF. At a particular material point X with an infinitesimal area dA with a unit normal vector N , we define the resulting trac- tion vector T in terms of the first Piola-Kirchoff (PK1) stress tensor P , T = PN. (8) In homogenization, the constitutive response on the struc- tural (or macroscopic) scale of a body depends strongly on the underlying microstructure. In the following, we consider 123 194 Computational Mechanics (2023) 71:191–212 a periodic microstructure � = [ − l1 2 , l1 2 ] × [ − l2 2 , l2 2 ] × [ − l3 2 , l3 2 ] . (9) Separation of length scales is assumed, i.e., established homogenizationprinciples applywithout the need for advanced modeling techniques, e.g., based on filtering [20] or higher- order continuum theories [21]. The objective of homogenization is the identification of the effective, macroscopic constitutive response P = 〈P〉� (10) of micro-heterogeneous materials for given F = 〈F〉�. Fur- ther, phase-averaged stresses, stress statistics, and the algo- rithmic tangent operator can be of relevance. The effective quantities must satisfy the Hill-Mandel macro-homogeneity condition P · F = 〈P〉� · 〈F〉� = 〈P(X) · F(X)〉� . (11) Periodic fluctuation boundary conditions of the form U(X) = (F − I)X + Ũ(X) (12) satisfy the Hill-Mandel condition (see, e.g., [22]) and have proven versatile and efficient. The tractions T (X±) are then anti-periodic for point-pairs X± ∈ ∂Ω± on opposing faces of the RVE �. The related function space for Ũ is referred to asV# which is a subset of the Sobolev space of weakly differentiable functions H1(Ω): V# = { Ũ ∈ H1(Ω) : Ũ(X+) = Ũ(X−), Ũ∗(X) = Ũ(mod(X + X∗,Ω)) ∈ H1(Ω) ∀X∗ ∈ Ω) } , (13) where mod is the spatial modulo operator. The homogeniza- tion problem (HOM) on � for prescribed deformation F reads: find Ũ ∈ V# (14) s.t.: Div (P) = 0 in � . (15) 3 FFT-based homogenization The solution of (HOM) (14)-(15) can either be obtained (in seldom cases) using analytical solution or by using discrete numerical techniques. The most widely used methods for computational homogenization in solid mechanics are cer- tainly the finite element method (FEM) [e.g., FE2, [23]] and FFT-based homogenization. The latter was proposed byMoulinec and Suquet [6] in the early 1990s. It is based on the Lippmann-Schwinger equa- tion in elasticity [24,25] and avoids both time-consuming meshing needed by conforming finite element discretiza- tions as well as the assembly of the related linear system. Therefore, thememory needed for solving the problem is sig- nificantly reduced compared with other methods. By virtue of the seminal Fast Fourier Transform algorithm [FFT, [26]], the compute time scales with O(n log(n)) where n is the number of unknowns, i.e., it is just slightly superlinear. For finite strain problems, the Lippmann-Schwinger equa- tion reads F = F − �0 : (P − C 0 : F), (16) with the Greens’ operator �0(·) = Grad ( G0Div (·) ) , (17) which relies on the solution operator G0 of a linear reference problem. Lahellec, Moulinec, and Suquet [27] proposed to solve (16) by the Newton-Raphson method and the linear Moulinec-Suquet fixed point solver. In contrast, Eisenlohr et al. [28] suggested using the Moulinec-Suquet fixed point iteration on the nonlinear Lippmann-Schwinger equation for finite strains directly. Kabel et al. [29] carried over the idea of Vinogradov and Milton [30], and of Gélébart and Mondon- Cancel [31] 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 con- jugate gradient (CG) method introduced by Schneider [32] for small strains. Further, the FFT-accelerated solution of reg- ular tri-linear hexahedral elements [FFT-Q1 Hex, see [12]] and in the closely related Fourier-Accelerated Nodal Solver [11] will be used in the sequel, for which the nonlinear CG method has a perfectly natural interpretation with the fun- damental solution acting as preconditioner. An overview of alternative solution methods can be found in [33]. Regarding the discretization, wewill compare the original discretization by Fourier polynomials of Moulinec-Suquet [34] with the HEX8R discretization of Willot [9], and the fully integrated HEX8 elements [12,35]. In order to consistently evaluate material laws for the staggered grid discretization [10,36], we apply the idea of the double fine material grid (DFMG) to large strains. In contrast to the regime for small strains, which can be efficiently imple- mented for isotropically linear elastic materials by a suitable averaging of the shear moduli, the DFMG for large deforma- tions requires multiple evaluations of the material routines, 123 Computational Mechanics (2023) 71:191–212 195 regardless of the complexity of the geometrically nonlinear material law. The technical details of the implementation can be found in the Appendix A. 4 Composites voxels One of the main advantages but—at the same time—a major disadvantage of FFT-based techniques is the constraint to reg- ularCartesiangrids:While they enable direct use of 3D image data, they are unable to capture the material interface accu- rately compared to interface-conforming discretizations. This is primarily due to thebinarizednature of themicrostruc- ture, which leads to the so-called staircase approximation of the interface. In order to capture the microscale effects suf- ficiently, a high grid resolution is hence needed, which–in turn–calls for high computational cost. In order to limit or even eliminate the staircase phenomenonwithout the need for a (distinct) grid refinement, so-called composite voxels were previously suggested in [15–17,37]. They enhance the exist- ing binary discretization using special voxels with effective material properties that depend on the phase volume fractions and the normal orientation N of the laminate. Consider such a composite voxel �e comprised of two phases denoted by�e± � �e, where+ and− stand for inclu- sion andmatrix 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 W±. Within �e, the material interface is approximated by a planeS e leading to a rank-1 laminate defined by the interface normal vector N ∈ R d in mate- rial configuration. The individual phase volume fractions are given by c+ and c− where c+ + c− = 1. In previous works, different rules of mixture have been suggested, namely the Voigt (•V) and Reuss (•R)1 estimates correspond to the upper and lower bounds of the mechanical response. Inspired by laminate theory, [38] suggests the definition of the effective elasticity tensor of the laminate C ‖ � implicitly through ( P + λ ( C ‖ � − λI )−1 )−1 = 〈( P + λ (C − λI)−1 )−1 〉 �e . (18) 1 In the nonlinear kinematic regime, these correspond to the Taylor and Sachs approximation, respectively. Here, the fourth order tensor P depends on the normal N via Pi jkl = 1 2 ( Niδ jk Nl + Niδ jl Nk + N jδik Nl + N jδil Nk ) − Ni N j NkNl (19) for i, j, k, l ∈ {1, . . . , d}. Note that in laminate theory, the states in the ± phase are piece-wise constant. Hence, the averaging translates into 〈(·)〉�e = c+ (·)+ + c− (·)−. (20) The laminate mixing rule considers the interface orienta- tion and introduces a parameter λ > 0 that must be larger than the leading eigenvalue of the individual stiffness ten- sors C±. Solving (18) for linear elastic materials involves 6 inversions of 6 × 6 matrices2 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. 4.1 Hadamard jump conditions for finite strain kinematics Following the established laminate theory, the following assumptions hold: • The deformation gradient in �e+ and �e− are related by a rank-1 jump along the normal orientation, i.e., for ã ∈ R d , �F�S e = F+ − F− = ã ⊗ N . (21) • The traction vector is continuous across the interfaceS e, i.e., �T�S e = �P�S eN = 0, T+ = T− . (22) The kinematic compatibility of the interface is character- ized by having a continuous deformation gradient tangential to the interface. This ensures that material point pairs remain identical. Condition (22) ensures that the interface is in static equilibrium. In order to enforce (21), we chose the following param- eterization3 of the deformation gradient tensors of the two 2 Symmetry is exploited. 3 (slightly different to that used, e.g., in [17]) 123 196 Computational Mechanics (2023) 71:191–212 material phases F± = F� ± 1 c± (a ⊗ N), (23) where a is related to ã in (21) by a scaling constant and F� is the prescribed deformation gradient on the composite voxel �e . The parameterization (23) preserves the volume average F� of the deformation gradient F over the compos- ite voxel �e according to 〈F〉�e = F� + ( c+ c+ − c− c− ) a ⊗ N = F� . (24) In the finite strain context, not only the average of F must be preserved, but also the total volume of the deformed material must remain constant under the rank-one perturbation: Lemma 1 The kinematic rank-1 perturbation on thematerial interface S e is volume preserving by construction. Proof Suppose A is an invertible n × n matrix and u, v are column vectors of size n. Then thematrix determinant lemma states that det ( A + uvT ) = (1 + vTA−1u) det (A) . (25) The deformation gradient F� is regular by definition since material inversion is not allowed. Thus, by setting A ← F� and u ← ±a/c±, v ← N the material Jacobian J± in the two phases can be computed using the matrix determinant lemma J± = det (F±) = ( 1 ± 1 c± aTF−T � N ) det (F�) . (26) Averaging J over the composite voxel �e, we obtain 〈J 〉�e = c+ J+ + c− J− = det (F�) = J� . (27) � To obtain the gradient jump vector a, the equilibrium of the tractions on the interface (22) must be granted; see also [17]: f (a) = T+ − T− != 0. (28) Examining the relative volume in �e±, additional constraints on J± emerge: J+ ! > 0 , J− ! > 0 . (29) These constraints imposed on (28) are essential to gain robustness, see Sect. 4.2. 4.1.1 Algorithmic implementation In the case of large strain kinematics, the system (28) is always nonlinear with no explicit solution available, in gen- eral. Hence, it must be solved iteratively, e.g., by using a Newton–Raphson scheme with the Hessian � f = d f (a) da = ( ∂ P+ ∂F+ ∂F+ ∂a − ∂ P− ∂F− ∂F− ∂a ) N. (30) The Hessian in index notation can be written in terms of the 2-tensors P± and F± as ( � f ) ik = NJ ( 1 c+ A i JkL+ + 1 c− A i JkL− ) NL , (31) where, A± = ∂ P± ∂F± = ∂2W± ∂F± ⊗ ∂F± . The constitutive tangent operatorA± can directly be obtained from the hyper- elastic potentials W±. Using a vector notation for F, P and the related matrix representation A± of A±, the Hessian gets � f = DT ( A+ c+ + A− c− ) D (32) with the matrix D defined via (a ⊗ N) → D a = ⎡ ⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ N1 0 0 0 N2 0 0 0 N3 N1 0 0 0 N2 0 0 0 N3 N1 0 0 0 N2 0 0 0 N3 ⎤ ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦ [ a1 a2 a3 ] . (33) Note the straightforward symmetric structure of thematrix � f (32). Within iteration [k], the Newton-Raphson update reads a[k] = a[k−1] − ( �[k−1] f )−1 f (a[k−1]) . (34) The full algorithm for the naive Newton–Raphson (NR) update is given in Algorithm 1. Note that each iteration requires a single d × d matrix inversion. Once convergence is achieved, i.e., once the tractions on the interface are in balance, the first Piola-Kirchoff stress is P� = c+P+ + c−P− . (35) 123 Computational Mechanics (2023) 71:191–212 197 The unconditionally symmetric effective stiffness matrix can be computed from δA = A+ − A− (36) ∂P� ∂F� = A� = ( c+A+ + c−A+ ) − δA D �−1 f DTδA . (37) Algorithm 1 Newton-Raphson (NR) algorithm solving for gradient jumps in composite boxels Require: F� Ensure: T+ = T− a ← 0 (or any admissible choice) � initialization F± ← F� ± 1 c± a ⊗ N P±, A± ← material model (F±) f ← (P+ − P−)N � compute initial residual E� ← | f | � initial error while E� > ε do � check convergence a ← a − �−1 f f � naive NR-update F± ← F� ± 1 c± a ⊗ N � update F± P±, A± ← material model (F±) f ← (P+ − P−)N � update residual E� ← | f | � update error end while 4.2 Robust algorithm via selective back-projection Unfortunately, a naive update of the gradient jump vector a ← a−�−1 f f sometimes leads to physically unacceptable iterates: the local volume within �± can get negative which is inadmissible, i.e., contradicting the constraints (29). Robustness can be enhanced by restricting the material Jacobian J± in both phases to be strictly positive according to (29). Previously, this problem was addressed in [16] by limiting the overall step width of the NR iterations if the constraint (29) is not met. We choose a different approach that can improve the con- vergence behavior by building on the matrix determinant lemma (25). It allows for rewriting J± as J± = ( 1 ± 1 c± aT ·F−T � N ) J� > 0. (38) We define β± c = ∓ c± ‖F−T � N‖ , mβ = F−T � N ‖F−T � N‖ (39) and, with help of mβ the orthogonal projectors M‖ β = mβ ⊗ mβ, M⊥ β = I − M‖ β . (40) The matrix determinant lemma implies that every NR iterate a must satisfy β+ c < aTmβ < β− c , (41) i.e., a condition that is trivial to check. Most importantly,mβ is independent of the current iteration, i.e., it can be precom- puted. If starting from a previously admissible point a0 = a[k−1] an inadmissible iterate a1 ← a0 + δa is detected, we com- pute βi = ai · mβ (i ∈ {0, 1}). Obviously β0 is admissible while β1 is not. We suggest to compute a selectively back- projected admissible coordinate β∗ as the mid point between the admissible β0 and the critical value βc ∈ {β+, β−} (the value which is exceeded by β1 is taken, see Algorithm 2): β∗ = β0 + βc 2 . (42) Other than the approach suggested in [16] our selective back- projection yields a different (admissible) iterate by adjusting only the part of a contributing to J±: a∗ = a0 + β∗ mβ = a1 + (β∗ − β1) mβ. (43) It leaves the part of a not interacting with the volume con- straint – that is, the part of a1 that is orthogonal to mβ – unaltered. Only the part of a1 that is co-linear with mβ is modified such that a valid iterate is obtained. Using the pro- jectors (40), this is equivalent to: M⊥ β a∗ = M⊥ β a1, M‖ β a∗ = β∗ mβ . (44) This induces that a large part of the NR update is preserved. No additional hyperparameters, line searches, and constitu- tive evaluations are required. Algorithm 2 Selective back-projection of inadmissible a1 for previously admissible a0 Require: inadmissible a1; a0; mβ, β± c Ensure: β+ c < aT∗mβ < β− c β1 ← aT1mβ , β0 ← aT0mβ � initialize if β1 ≤ β+ c then � get critical value βc βc = β+ c else if β1 ≥ β− c then βc = β− c end if β∗ ← (β0 + βc)/2 � get admissible β value a∗ ← M⊥ β a1 + β∗ mβ � update a accordingly 123 198 Computational Mechanics (2023) 71:191–212 The increased robustness due to the selective back- projection is of utmost relevance for actual finite strain simulations, particularly in FE2 like settings. This holds a fortiori if composite voxels with rather small phase volume fractions go along with pronounced phase contrast and high load increments, see Sect. 9. As mentioned earlier, our algorithm shares similarities with the backtracking algorithm proposed in [16]. Both approaches use the matrix determinant lemma to check for inadmissible iterates. The backtracking of [16] is employed as a line search strategy to compute an optimal step size in the admissible range, leading to scalar backtracking of a and a damped NR method. This induces additional line search parameters as well as additional function evaluations. Con- trary to that, our selective back-projection algorithm requires no extra constitutive evaluations and ensures an admissible update unconditionally in the absence of additional hyperpa- rameters. 5 Identification of the normal vector for composite voxels and boxels 5.1 Existing procedures and complications Composite voxels and boxels (non-equiaxed voxels) require the normal vector N alongside the phase volume fractions c±. While c+ and c− are easy to compute from averaging over the composite voxel or boxel, identifying the normal vector N is not without issues. It is also evident that the quality of the normal information is critical to the performance of the com- posite boxels. The original proposal of Kabel and colleagues [15] uses the direction of the line connecting the barycen- ter of the dominant of the two phases within the composite voxel/boxel against the center of the voxel as an approxima- tion of the normal vector, see Fig. 1. A trivial computation shows that this is equivalent to having the normal defined via the connecting line of the barycenters: Ñ = c− − c+ |c− − c+| . (45) In the present work, the normal is assumed to point out of the inclusion phase (indexed by +). The approach (45) has several advantages to it: The imple- mentation is rather trivial, the computation is rapid, and it is easy to guarantee that the normal N points from phase 1 (denoted inclusion phase in the sequel) into phase 0 (denoted matrix phase in the sequel). However, the use of the barycen- ters c± for the normal detection is not without issues, as can be seen from the examples shown in Fig. 1. At low volume fractions, the normal direction will depend only on the posi- tion of the phase rather than on the actual interface (compare (c) and (d) in Fig. 1), and the method does not work for boxels that differ from equiaxed composite voxels by hav- ing (sometimes pronounced) aspect ratios, see Fig. 2. This occurs if the number of voxels inside the composite boxel varies along the different edges (with fine-scale voxels being equiaxed, e.g., Fig. 2, top). For instance, this can be useful in order to reduce the resolution in pseudo-unidirectional fiber reinforced materials; see Sect. 6.5 for examples. In order to allow for an improved normal computation over (45), the authors suggest a novel strategy that can be broken down into two steps: N.1 compute an indicator for the interface between the phases using a discrete Laplacian resulting in discrete weights on the fine grid; N.2 within composite boxels that contain a material interface define the normal vector via a minimization problem. The individual steps are described below. A free python implementation is available via an open access software accessible via GitHub ( [39]), including the option to pro- cess HDF5 files easily [40] (see also the documentation and tutorial given in the repository’s Jupyter notebook). 5.2 Interface indication via a discrete Laplacian In the following, the Laplace filter commonly used in image processing is introduced. It can be employed for edge detec- tion. For simplicity, the Laplace stencil is defined by building on a first-order finite difference gradient operator along each coordinate axis. Therefore, the 3D stencil illustrated in Fig. 3 is utilized. In 1D for grid spacing h > 0, the derivative at position xi and the second derivative (gathered from the subsequent application of the first derivative) read ∂ f (xi−1) ∂x ≈ f ′ i−1 = fi − fi−1 h , (46) ∂ f (xi ) ∂x ≈ f ′ i = fi+1 − fi h , (47) ∂ f 2(xi ) ∂x2 ≈ 1 h ( f ′ i − f ′ i−1) = 1 h2 ( fi−1 − 2 fi + fi+1) . (48) Setting f j = f (x j ) to value 1, if the inclusion phase is found at x j and to 0 otherwise, the Laplacian will be 0, if and only if all pixel values within the stencil are identical. Therefore, only point triples hostingmore than one phasewill lead to a nonzero Laplacian. Further, pixels with positive and negative valueswill be found,which represent inside and out- side voxels related to the interface, respectively. For 2D and 3D boxels, the stencil is composed of uniaxial Laplace sten- 123 Computational Mechanics (2023) 71:191–212 199 (a) (b) (c) (d) Fig. 1 Normal direction following [15]: Ñ from (45) is the direction between the barycenters of the phases for different examples cils along the coordinate axes, using the respective h value corresponding to the fine-scale boxel dimension along the respective direction (hx , hy, hz), see also Fig. 3. The Laplace stencil will be weighted by the boxel volume. Thereby, the scalar factor will get a neat physical interpretation: ri = r(li ) = l j lk li (i �= j, i �= k, j �= k). (49) The ratio ri expresses the ratio of the area of the boxel face with normal direction ei divided by the boxel dimension along direction ei . After application of the stencil S to the image, the absolute value is taken to gain the discrete weights wi jk = |S ∗ χ |i jk (i ∈ {1, n1}, j ∈ {1, n2}, k ∈ {1, n3}) . (50) The convolution is carried out in the Fourier domain, which automatically enforces periodic boundary conditions for the normal detection of inclusions crossing the edges of the computational domain. Of course, also a direct computation would be feasible at a comparable compute cost. Theweights are now processed for each composite boxel attributed with multiplematerial phases, i.e., with volume fraction 0 < c− < 1. On each of these boxels, a weighted least squares problem is set up in order to identify the normal orientation. Therefore, on the boxel �B the second moment tensor M = ∑ {i, j,k}∈�B wi jkxi jk ⊗ xi jk ∈ Sym+(R3) (51) is computed, where xi jk denotes the coordinates of the voxel with discrete coordinates (i, j, k)4. Next, the eigenvector matching the smallest eigenvalue of M is used as the ini- tial normal N . In a second step, the direction of the normal N is identified. Therefore, the vector connecting the barycen- ters of the material phase within the boxels Ñ given by (45) 4 It is assumed that xi jk is relative to the barycenter of ��. Fig. 2 Normal direction Ñ from (45) following [15] for a 2D boxel (16×6 fine-scale grid, equiaxed fine scale); the true normal N0 is pro- vided for comparison of the original approach [15] is considered to identify the proper sign of the proposed N : N ← { N if Ñ · N > 0 −N else. (52) Thereby, the normal is guaranteed to point out of the+ phase. The example shown earlier in Fig. 2 has been used with our approach, see Fig. 4. The interface voxels are shown in red with darker tones denoting increased weights. Note the good correlation of the analytically defined normal N0 used to generate the discrete image and the normal reconstructed using the proposed procedure. The collinearity of the two vectors is N0 · N = 0.99998, where the interface detection was effected using matching padding for the analytical nor- mal. This compares against N0 · Ñ = 0.7761 for the normal shown in Fig. 2 obtained from (45). Another comparison is shown in Fig. 5. Here an input image consisting of 2563 voxels containing a centered sphere of radius r = 0.4 L (with L denoting the edge length of the cube) is coarsened using composite voxels of size 323, i.e., the information is compressed by a factor 323=32768 yielding a coarse scale resolution of just 83. For each com- posite voxel the normal is computed using either the approach from (45) [15] (Fig. 5a) and the approach presented in this Section (Fig. 5b). In this visualization, the actual facets are reconstructed from the normal and the volume fraction c+. 123 200 Computational Mechanics (2023) 71:191–212 Fig. 3 3D stencil used for the Laplacian (here h1 = h2 = h3 = 1 for simplicity) Fig. 4 Normal vector for the same example as in Fig. 2 but using the proposed algorithm; interface voxels are red; cS is the barycenter of the interface By the metric of vision, the normals of the old approach graphs are less regular. Obviously, the facets reconstructed from these are also not leading to an accurate approximation of the curved surface of the sphere (e.g., gaps/overlaps). Our procedure leads to visually more accurate normals. This is confirmed by the planar facet reconstruction that is almost entirely devoid of gaps and overlaps. Next, we have repeated the previous comparison by scaling down from 2563 to 8×16×32, i.e., introducing non- equiaxed composite boxels of size 32×16×8. The top views of the concurrent approaches for the normal computation are compared in Fig. 5. The difference when using boxels is massive, with large gaps showing in Fig. 5c for the old method,while basically nogaps andminimal overlap is found in Fig. 5d for the improved normal identification. The results from Figs. 4 and 5 are consistent. They emphasize the rele- vance of using a dedicated scheme for normal identification. This is evermore so true in the presence of small volume fractions and/or non-equiaxed boxels. 6 Numerical examples 6.1 Material models and loading conditions In this section, we will focus on the key aspects and novelties of the proposed composite boxel approach for mechanical applications at finite strains. We assume that the materials obey a compressible Neo-Hookean material model W = 1 2 λ(ln J )2 − μ ln J + 1 2 μ(tr(C) − 3). (53) The Young’s modulus E± and the Poisson ratios ν± for the inclusion and the matrix phase are denoted by the respec- tive subscripts. We have deliberately chosen a rather distinct contrast by setting the synthetic parameters E+ = 10 GPa, ν+ = 0.3 E− = 1 GPa, ν− = 0.0. (54) Note that Young’s modulus has a distinct contrast of 10, exceeding that of literally all practical metal matrix com- posites, while the difference in the Poisson ratio is also pronounced5. Elevated contrast in the Poisson ratio is usually detrimental to the convergence of FFT-based schemes as it alters the collinearity of the stiffness of the contained phases, see [11] for a convergence study. The macroscopic deformation gradient F imposed on the RVE in the homogenization problem (HOM) is chosen to be 50% pure shear in component Fxy F = ⎛ ⎝1 1/2 0 0 1 0 0 0 1 ⎞ ⎠ . (55) The authors also emphasize the arbitrary selection of the materialmodels (trying to trigger elevatedmaterial contrasts) and of the imposed kinematic loading. In view of repro- ducibility, we have opted for a simple unambiguous model and a substantial loading consideringmultiscale applications. 6.2 Spherical inclusion Note that all of the following examples are truly 3D, although some slice views could lead to the misinterpretation of a 2D problem. First, we look at a benchmark problem with a spherical inclusion of radius R = 0.4 embedded in a 3D matrix. We aim to identify the impact of a more accurate identification (Sect. 5) for both equiaxed and non-equiaxed composite box- els. The influence on the phase averages is also studied. We 5 This is equivalent to a contrast of 25 in the bulk modulus K and a contrast of 7.69 in the shear modulus G. 123 Computational Mechanics (2023) 71:191–212 201 (a) (b) (c) (d) Fig. 5 (Top view) Comparison of the established normal computation [15] a, c and the proposed approach b, d with composite boxels of size 32 × 32 × 32 and 32 × 16 × 8 respectively. Input image of size 2563 start from a fine-scale microstructure of 2563 voxels which is down-scaled to a resolution of 323 voxels (downscale factor 512). The coarsened discretization comprises 2624 (∼ 8%) composite voxels. In each composite voxel, the normal ori- entation N and the local phase volume fractions c± are computed via Sect. 5 and using the established method of [15]. The homogenization problem is then solved via Fourier- Accelerated Nodal Solvers (FANS, [11]) using the 8-noded hexahedral finite elements. The convergence criterion is a relative tolerance of 10-10 with respect to the l∞-norm of the nodal force residual vector. In Fig. 6, the converged solution is shown at the center slice Z = 0. We employ a post-processing procedure that gathers individual field data for each phase and visualizes it using the planar interface with orientation N . It is obvious that this leads to a much smoother solution close to the interface, even within the individual phases. To our knowledge, a comparable visualization has not been used in previous studies. Further, we have visualized the tractions in the deformed configuration (Fig. 6c), which are unavail- able in both non-conforming discretizations as well as in usual conforming FE discretizations. Regarding the accuracy of the effective stress tensor P and the phase averages P±, a study was performed which compares the 2563 reference solution using the procedure suggested in [9] for five different ComBo discretizations using reduced integration CG-FANS [11] (FANS HEX8R) and using the normals from (45) matching previous studies as well as the improved normals from Sect. 5. The results are summarized in Table 1. Note that the average stress only reflects part of the actual accuracy of the solver as it neglects local field fluctuations, which are examined in the following examples, see Sect. 6.3. It can be noted that the errors in P of all composite discretizations are below 1.15% for the nor- mals of [15] and 0.78% for the normals suggested in Sect. 5, respectively, i.e., the improved normals reduce the error in Table 1 Comparing averaged 1st Piola-Kirchoff stresses in the spher- ical inclusion problem (Sect. 6.2): different resolutions, normals obtained via [15] and Sect. 5, with errors against a reference solution (2563) based on [9]; ComBo solutions are obtained using FANS with reduced integration [11]; bold highlights better result (old normals vs. proposed normals); Error defined as the relative Frobenius norm w.r.t reference solution Resolution Downscale Normals via [15] Normals via Sect. 5 Error (%) Error (%) P P− P+ P P− P+ 8 × 8 × 8 32768 1.145 0.263 3.537 0.771 0.177 2.361 16 × 16 × 16 4096 0.160 0.042 0.501 0.053 0.028 0.172 8 × 16 × 32 4096 1.586 0.282 4.609 0.177 0.064 0.584 32 × 32 × 32 512 0.072 0.019 0.219 0.070 0.018 0.212 16 × 32 × 64 512 0.801 0.156 2.335 0.100 0.019 0.262 all averaged stresses by approximately 30% for equiaxed composite voxels and up to 1000% for the non-equiaxed composite boxels. Remarkably, these improvements come at no additional computational expense, but they owe only to the more accurate orientation information. 6.3 Composite boxels and local solution field quality The straightforward implementation of composite boxels into various FFT-based schemes makes them truly versatile. Here, we present some of the most popular methods used in tandem with composite boxels. Most importantly, the com- posite boxel method can be used to replace any call to a constitutivemodel, independent of the discretizationmethod. In Fig. 7, we consider a random polyhedral inclusion sur- rounded by thematrixmaterial. Thefine-scalemicrostructure is generated at a resolution of 2563. The ComBo discretiza- tion is 323 with normals gained from Sect. 5, yielding a 512 times smaller problem. We further go on to compare the full 123 202 Computational Mechanics (2023) 71:191–212 (a) (b) (c) Fig. 6 Composite Boxel (ComBo) solution using CG-FANS [11] with full integration for the spherical inclusion example of Sect. 6.2: spatial configuration at slice Z = 0 field solutions for the EXY component of theGreen-Lagrange strain at the center slice (Z = 0) for various popular FFT- based solution strategies. A reference solution is computed on the original fine-scale problem (without any composite boxels). The solution using FANS with full integration (8 Gauss points per element) finite elements (HEX8) is shown next to the reference. It exhibits no staircasing and no artifacts even close to the interface. The solution quantitatively captures the behavior well in the interior of the phases and near the phase boundaries compared to the reference solution, even for this coarse discretization. Stereotypical of the fully inte- grated HEX8 elements, the response is a tad stiffer, which is reflected in the phase-wise averages, Table 2. The solution obtained by the staggered grid approach, despite its advan- tages in an algorithmic sense, lacks symmetry in the physical location of the strain components in each voxel. Hence, only an interpolated measure can be obtained at the boxel cen- ter for visualization. This ambiguity, unfortunately, leads to poor results close to the material interface despite the use of the ComBo discretization, which is also reflected in the homogenized quantities. This issue can partially be over- come by the proposed double-fine material grid (DFMG, see AppendixA) approach at the cost of extramaterial law evalu- ations compared to the original staggered grid approach. The DFMG approach also interpolates the field quantities to the boxel center, which causes blurring/smoothing everywhere, including the material boundaries: local solution field accu- racy is sacrificed, althoughDFMGoutperforms the staggered grid approach, and it preserves existing symmetries. FANS HEX8R, i.e., with reduced integration (1 Gauss point per ele- ment), is numerically similar to the HEX8R discretization by Willot [9], as shown by [10]. The FANS HEX8R solution suffers from hourglassing, but the amplitude of the hourglass modes is greatly diminished in the presence of composite Table 2 Comparing averaged 1st Piola-Kirchoff stresses in the Poly- hedron problem (323) (Sect. 6.3): different FFT based methods: with errors against a reference solution (2563) based on [9]; Error defined as the relative Frobenius norm w.r.t reference solution Error (%) Method P P− P+ Staggered Grid 1.219 0.202 16.993 DFMG 0.661 0.103 9.147 FANS HEX8 0.142 0.024 1.947 Willot 0.065 0.011 0.903 FANS HEX8R 0.049 0.008 0.675 Moulinec-Suquet 0.026 0.005 0.349 boxels compared to when no composite boxels are in use. The convergence behavior of FANSHEX8R is very much on parwith FANSwith fully integratedHEX8 elements. Finally, we also employ the originalMoulinec-Suquet scheme,which has been extensively studied in the literature. It suffers from spurious oscillations, although predicting the homogenized quantities very well. 6.4 Selective back-projection in practice The selective back-projection introduced in Sect. 4.2 is inves- tigated for the polyhedral microstructure of the previous Sect. 6.3. First, we demonstrate the influence it has on the convergence of the NR scheme inside of a single critical voxel, see Fig. 8. The traction balance residual vs. material Jacobian of phase − is plotted in a particular case when back-projection is employed. In the case discussed, the volume fraction of the inclusion phase + is 99.625%. The Newton-Raphson algo- rithm 1 starts with an admissible initial guess (at k = 0), but the naive NR-update would push the iterate a[1] to the 123 Computational Mechanics (2023) 71:191–212 203 Fig. 7 Solutions for the Green-Lagrange strain EXY obtained by different FFT-based schemes for the polyhedral inclusion problem (Resolution: 323) and comparison against reference solution (Resolution: 5123) at slice Z = 0 Fig. 8 Convergence behavior of the Newton–Raphson algorithm with back-projection inadmissible domain, which is highlighted in red in Fig. 8. The solid blue line tracks the continuous intermediate resid- ual between the current iterate a[k] and the subsequent iterate a[k+1] along a line search parameter. Note that if either one of J± tends to 0, the residual tends to ∞, which is character- istic of log based hyperelastic strain energies [41]. It is also noteworthy that, although the solution lies in between a[0] and a[1] in this plot, the continuous residual does not go to zero anywhere. This is due to the projection onto the J− axis (corresponding to the β axis modulo scaling) which cannot account for incorrect components M⊥ β a [k] – i.e., a mere line search cannot suffice, in general. Using the selective back- projection algorithm 2 yields a valid iterate a[1]. From there on, the NR algorithm 1 converges to a physically feasible solution within six iterations up to machine precision, and in practice, one could stop after just four iterations. Each iter- ation of our NR algorithm equates to a single evaluation of the constitutive law for either phase, independent of whether selective back-projection is needed. Despite striking similarities with the algorithm proposed in [16], we developed our scheme independently since we observed that it was needed during the simulations.While the authors of [16] state that “The projection and backtracking steps occur only rarely”, we would like to emphasize that a single inadmissible iterate can break the entire simulation (e.g. leading to negative J± used in log-energies). Further, we have investigated the percentage of the composite boxels that require selective back-projection as a function of the loading for the polyhedral inclusion, see Fig. 9. 123 204 Computational Mechanics (2023) 71:191–212 Fig. 9 Percentage of boxels (Problem c.f. Sect. 6.3) with failed naive N-R updates for mixed displacement gradient loading conditions It turns out that our algorithm is needed in a majority of the composite voxels under certain loading conditions. This becomes critical in actual two-scale simulations, where spuriously high loadings are frequently observed in isolated points, especially in the compression regime. Ona side note,wewould like to emphasize that in our tests, the old normal direction cf. [15] further increased the number of failed iterates. Hence, using the normal detection from Sect. 5 alongside Algorithm 2 improves on the robustness in both regards. 6.5 Short fiber reinforcedmicrostructures In this Section, we investigate fiber-reinforced composites with global fiber directional affinity to show the effective- ness of non-equiaxed composite boxels in certain use cases. We consider the twomicrostructures shown in Fig. 10a and b referred to as [FPR-1] and [FRP-2] respectively. Microstruc- ture [FRP-1] is a fibrous microstructure with almost aligned fiberswith afiber volume fraction of∼ 4%primarily oriented along the X -axis. This problem is down-scaled to varied equiaxed and non-equiaxed resolutions starting from a fine- scale image comprising 5043 voxels tabulated in Table 3. Similarly, themicrostructure [FRP-2] hosts 150 fibers almost aligned along the X -axis with a fiber volume fraction of ∼ 15%. The fibers have a length of 120μm, a diameter of 12μm, and a minimum fiber distance is about 2μm. This problem is down-scaled to varied equiaxed and non-equiaxed resolutions starting from a fine-scale image comprising 2403 voxels tabulated in Table 4. The overall homogenized first Piola-Kirchoff stress and its phase-wise counterparts are compared towards a reference solution computed without the use of composite boxels and using the FFT solver proposed in [9]while the coarse-grained models are solved using FANS HEX8R. Figure 10c and d show that theComBodiscretization leads to reasonable local stress fields despite the use of massively anisotropic grids. Thereby, the resolution of the simulation can adapt to the aspect ratio of the fibers, allowing for dis- tinct computational gains while, at the same time, the lateral resolution can remain sufficiently fine to gain (a) accurate ori- entation data (Sect. 5) and (b) to separate individual fibers. Both can be achieved without growing the number of degrees of freedom. For instance, [FRP-1] is resolved 3.5 times coarser along the fiber axis, yielding a speed-up and mem- ory savings in the same order against an equiaxed grid. For [FRP-2] several ComBodiscretizations are compared against a reference solution regarding P, P± in Table 3. The overall accuracy was better than 1% with improvements indepen- dent of the boxel aspect ratio as the number of boxels grows. The accuracy within the inclusion phase is impressive, given that the fibers make up only ∼ 4% of the material. Simi- lar trends are observed in Table 4 for the [FRP-2] problem, where the use of non-equiaxed boxels becomes vital in the need to resolve very small fiber distances.Much coarser non- equiaxed resolutions still perform better than corresponding equiaxed resolutions outlining the need to resolve the lat- eral plane of the fibers properly. It can also be observed that modest improvements are seen when the resolution along the fibers is increased, while much more substantial improve- ments are observed when the fiber lateral plane is refined. We also extract the Von-Mises Cauchy stress: overall and phase-wise histograms for the [FRP-2] problem as shown in Fig. 11, in order to judge the quality of the full field solution. The stress distributions for the equiaxed 483 (Downscale : 125) and the non-equiaxed 30 × 60 × 60 (Downscale: 128) is compared against a reference solution 2403. Although by the metric of vision, we can observe a better agreement with the non-equiaxed resolution against the reference solution, the corresponding cumulative squared deviations of the dis- tributions from the reference stress histogram are plotted, which unequivocally states that the non-equiaxed resolution has lower error across the board and thus captures the over-all field solution better. 7 Résumé 7.1 Summary We present an extension of the composite voxel/boxel (anisotropic voxel) approach of [15] towards finite strain hyperelasticity for FFT-based homogenization schemes sim- ilar to [16]. The foundations of FFT-based homogenization are recalled in Sect. 3, and a detailed description of the doubly-fine material grid which can rule out some issues of the staggered grid approach [10] regarding the local field accuracy, is outlined in Appendix A. 123 Computational Mechanics (2023) 71:191–212 205 (a) (b) (c) (d) Fig. 10 Fiber-reinforced composites: [FRP-1] ∼ 4% fiber volume fraction (a) [FRP-2] ∼ 16% fiber volume fraction (b); simulation results using FANS HEX8R with ComBo discretization 18×63×63 for [FRP-1] and 30×60×60 for [FRP-2] are shown in (c) and (d) Table 3 Comparing averaged 1st Piola-Kirchoff stresses in the [FRP-1] problem (Sect. 6.5): different equiaxed and non-equiaxed ComBo res- olutions, with errors against a reference solution (5043) based on [9]; ComBo solutions are obtained using FANS with reduced integration [11]; Bold highlights better result (equiaxed vs. boxels): Error defined as the relative Frobenius norm w.r.t reference solution Equiaxed voxels Non-equiaxed boxels Res Error (%) Downscale Error (%) Res Aspect ratio P P− P+ P P− P+ 213 0.894 0.148 13.916 13824 12348 0.654 0.110 10.205 18 × 242 4 : 3 : 3 243 0.570 0.097 8.913 9261 9072 0.461 0.080 7.219 18 × 282 14 : 9 : 9 283 0.351 0.061 5.503 5832 5488 0.298 0.051 4.662 18 × 362 4 : 2 : 2 363 0.189 0.033 2.975 2744 3024 0.168 0.030 2.650 24 × 422 7 : 4 : 4 423 0.119 0.021 1.871 1728 1792 0.075 0.017 1.222 18 × 632 7 : 2 : 2 633 0.035 0.007 0.555 512 504 0.036 0.006 0.564 36 × 842 7 : 3 : 3 843 0.037 0.006 0.579 216 224 0.049 0.007 0.744 36 × 1262 7 : 2 : 2 123 206 Computational Mechanics (2023) 71:191–212 Table 4 Comparing averaged 1st Piola-Kirchoff stresses in the [FRP-2] problem (Sect. 6.5): different equiaxed and non-equiaxed ComBo resolutions, with errors against a reference solution (2403) based on [9]; ComBo solutions are obtained using FANS with reduced integration [11]; the italic highlights better results than the corresponding equiaxed voxels resolution: Error defined as the relative Frobenius norm w.r.t reference solution Equiaxed voxels Non-equiaxed boxels Res Error (%) Downscale Error (%) Res Aspect ratio P P− P+ P P− P+ 300 0.870 0.186 3.759 20 × 482 12 : 5 : 5 250 0.783 0.169 3.392 24 × 482 2 : 1 : 1 403 1.110 0.234 4.806 216 200 0.726 0.157 3.149 30 × 482 8 : 5 : 5 192 0.375 0.081 1.586 20 × 602 3 : 1 : 1 160 0.299 0.070 1.283 24 × 602 5 : 2 : 2 483 0.688 0.146 2.982 125 128 0.238 0.059 1.041 30 × 602 2 : 1 : 1 108 0.254 0.039 0.973 20 × 802 4 : 1 : 1 96 0.228 0.055 1.004 40 × 602 3 : 2 : 2 603 0.221 0.048 0.961 64 72 0.111 0.027 0.451 30 × 802 8 : 3 : 3 The detailed algorithmic treatment of the composite vox- els/boxels in Sect. 4 yields low-, i.e.,d-dimensional nonlinear equations to be solved with explicit Hessians being pro- vided for infinitesimal and finite strain problems; see also the cheat sheet in Appendix 5. The algorithmic tangent oper- ator of the composite voxels is provided too, and it has a sleek representationwith a simplistic implementation. A cru- cial ingredient in composite boxel finite strain simulations is the back-projection scheme described in Algorithm 2. It ensures the admissibility of the deformation in either of the laminate phases at negligible computational overhead but much-increased robustness. When examining composite voxels and boxels, some issues with the normal detection [15] were found. In Sect. 5, an improved algorithm for normal identification is suggested, which leads to considerable improvements in the lami- nate orientation within the composite voxels. The proposed algorithm can also process composite boxels (ComBo) char- acterized by non-equiaxed coarsening,whichwas impossible using the approach by [15]. The procedure was shown to yield accurate normals for different microstructures. An open-source python implementation with examples can be accessed through the GitHub repository [39]. It also features 3D tools for visualization and a tutorial demonstrating the usage. In Sect. 6 a variety of different microstructures are sim- ulated using different FFT-based solvers, different normal detection procedures, and using different coarse-grained res- olutions. The results demonstrate that the local fields using theComBodiscretization are closelymatching full resolution solutions. FANS HEX8 and HEX8R [11,12] were found to yield the smoothest representation of the local fields. Despite the tendency of HEX8 to overestimate the stresses, this dis- cretization has the smoothest stress fields. Moreover, the suggested normal detection was shown to provide notably improved accuracy on the overall stress response as well as for the phase-wise averaged stresses (see Tables 1, 2, 3 and 4). The improvement for actual boxels was even more notable. Further, the stress statistics for equiaxed and non- equiaxed resolutions with similar downscale factors hint at the improved quality of the local stresses for the same down- scale factor. Additionally, the normals were shown to deliver convergence of the effective stresses that depends mainly on the number of DOF, i.e., the overall amount of coarse- graining, see Tables 3 and 4. Computational savings of 2000 and beyond at errors around 1% in the phase-averaged stresses are observed. 7.2 Discussion First up, the authors are thoroughly convinced that composite boxels have proven to be a valuable addition to many estab- lished FFT-based homogenization schemes. They allow for impressive computational savings in CPU time and memory (factor 2000 and beyond) at a modest—if any—sacrifice in accuracy. The accuracy was further improved by using the novel strategy for the normal detection from Sect. 5. It leads to a reduction of approximately 30% in the relative errors of the effective stress P and its phase-wise counterparts P± even for relatively smooth and simple microstructures. In the case of anisotropic boxels, more distinct improvements in the error were found. Surprisingly, in the presence of com- posite boxels, the number of DOF of the system seems to be the primary influence factor regarding the accuracy of the simulation even when pronounced boxel anisotropy is considered. A key advantage of non-equiaxed boxels is that pseudo-unidirectional fiber separation can be granted with- out growing the number of degrees of freedomof the problem while retaining accuracy. We think that this can leverage the simulations in cer- tain fields, e.g., for discontinuous short fiber composites with pronounced aspect ratios (e.g., 20 and beyond) and pseudo unidirectional fiber orientation. We are also convinced that making the source code for the normals freely available [39] 123 Computational Mechanics (2023) 71:191–212 207 Fig. 11 Normalized weighted distribution of the Von-Mises Cauchy stress (top) and the squared cumulative error (bottom) of the distribution w.r.t a reference (2403) solution could help in rendering composite boxels an attractive choice in academia and industry. In the future, we are confident that more refined comparisons of the actual solution fields of ComBo and high-resolution simulations will lead to further insights regarding the accuracy and overall efficiency. By the introduction of the doubly-fine material grid for the staggered grid discretization [10], considerable improve- ments with respect to the quality of the local fields were observed. However, these come at the expense of a distinct rise in the number of constitutive evaluations and little gain regarding the overall homogenized response. Therefore, this method is probably best suited when local solution fields are sought-after. In this regard, FANS [11], or the equiva- lent FFT-Q1 Hex [12] show the most confidence-inspiring results: local fields are smooth and match the reference solu- tion closely; hour-glassing is less distinct (HEX8R) or absent (HEX8) than in the almost identical discretization of [9]. It is important to state that the use of the ComBo discretization comes without issues. The application of composite voxels/boxels in finite strain homogenization problems revealed that particular care should be taken in view of considering physical constraints: The selective back-projection algorithm (Algo- rithm 2) demonstrates that robustness can be gained and (sometimes unrecognized) physically questionable iterates might occur in practice. We are confident that the presented algorithm is a leap forward. Despite this improvement, the realizable loading can still be limited, particularly when the phase volume fractions within the composite boxels tend towards 0 or 1 and the contrast in stiffness is pronounced.This can imply that the deformation gradient can approach critical states not just in the laminate but on the overall ComBo voxel (denoted F�). Finding further refinements to the simulation scheme could further boost robustness, giving rise to future research topics. An important message is also given by demonstrating the usefulness of ComBo discretizations for a rich set of differ- ent discretizations: FANS HEX8(R)/FFT-Q1 Hex [11,12], DFMG, staggered grid [10] the rotated grid scheme ofWillot [9], and the classicalMoulinec–Suquet scheme [6,7] were all used with success and building on the same implementation. The authors would like to emphasize that the approach is, however not limited to FFT-based schemes: regular Finite Element and Finite Difference schemes could use them too, yielding potential benefits without the intrusiveness of, e.g., the extended finite element method (X-FEM, e.g., [42]). Related to the recent progress reported by [43], the exten- sion of our framework for interface mechanics in the small and finite strain setting is a promising route. Major bene- fits due to the improved normal orientations from Sect. 5 are expected: Both, the local solution fields as well as the (thereby influence) interfacial tractions are assumed to gain in accuracy. Last, an equivalent to composite boxels for materials with more than two phases is urgently needed, e.g., in order to deal with polycrystals. 123 208 Computational Mechanics (2023) 71:191–212 Supplementary information The normal detection algorithm is available from [39]. Acknowledgements Funded by Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy - EXC 2075 - 390740016. Contributions by Felix Fritzen are funded by Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) within theHeisenberg program -DFG-FR2702/8 - 406068690. We acknowledge the support by the Stuttgart Center for Simulation Science (SimTech). Funding Open Access funding enabled and organized by Projekt DEAL. Declarations Conflict of interest The authors declare no potential conflict of interest. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adap- tation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indi- cate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, youwill need to obtain permission directly from the copy- right holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Appendix A A finite difference discretization on a staggered grid In the sequel, we describe the consistent formulation of the staggered grid discretization [10] for the geometrically non- linear case [36]. Similar to the small strain case, the diagonal and off-diagonal components of the deformation gradient are located at different positions; compare Fig. 12. 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 n1, n2, n3 and consider a regular peri- odic grid consisting of n = n1n2n3 cells, each a translate of [0, h1] × [0, h2] × [0, h3] with h j = l j n j ( j = 1, 2, 3). Let i+, j+ and k+ denote the location of the staggered grid coor- dinates, i.e.,a+ = a+ 1 2 fora = i, j, k. For the (i, j, k)-cell the coordinates of the displacement vector (U1, U2, U3) are located at the cell face centers, i.e., U1[i, j, k] is located at the coordinate (ih1, j+h2, k+h3), U2[i, j, k] lives at (i+h1, jh2, k+h3) and U3[i, j, k] is associated to (a) (b) Fig. 12 Placement of the strain and displacement variables in 2D: a Location of strain and displacement variables. b The variables’ grid (gray box) and the doubly fine material grid (orange box) (i+h1, j+h2, kh3). The situation is displayed in Fig. 12a, where we have restricted ourselves to the 2D case for clarity. The diagonal components FaA[i, j, k] (a, A ∈ {1, 2, 3}) of the deformation gradient F are positioned on the cell cen- ters (i+h1, j+h2, k+h3), whereas the off-diagonal strains F23, F32[i, j, k], F13, F31[i, j, k] and F12, F21[i, j, k], are locatedon the corresponding edgemidpoints (i+h1, jh2, kh3), (ih1, j+h2, kh3) and (ih1, j, (k+h3). For visualization, we again refer to Fig. 12a. Displacements and gradients are connected by central dif- ference formulae, where periodicity is understood implicitly. More precisely, introduce forward and backward difference operators on a scalar discrete field φ : Vn → R by the for- mulae D± j φ[I ] = ±φ[I ± e j ] − φ[I ] h j (A1) 123 http://creativecommons.org/licenses/by/4.0/ http://creativecommons.org/licenses/by/4.0/ Computational Mechanics (2023) 71:191–212 209 for I ∈ Vn = {0, 1, . . . , n1 − 1} × {0, 1, . . . , n2 − 1} × {0, 1, . . . , n3 − 1}. Then we introduce the gradient operator Grad (U) = ⎡ ⎣ D+ 1 U1 D− 2 U1 D− 3 U1 D− 1 U2 D+ 2 U2 D− 3 U2 D− 1 U3 D− 2 U3 D+ 3 U3 ⎤ ⎦ (A2) giving rise to the deformation gradient F = I + Grad (U) associated to a periodic displacement field U = (U1,U2,U3) : Vn → R 3, Similarly, there is a divergence operator, turning P : Vn → R 3×3 into Div P = ⎡ ⎣ D− 1 P11 + D+ 2 P21 + D+ 3 P31 D+ 1 P12 + D− 2 P22 + D+ 3 P32 D+ 1 P13 + D+ 2 P23 + D− 3 P33 ⎤ ⎦ . (A3) The stress variable PaB is located on the same position as the corresponding FaB . The constitutive law P(F) is defined piecewise on the regular voxel grid. For a typical material cell, shaded in gray in Fig. 12b, it is not clear how to apply the material law for off-diagonal strains, as the function P(F) may be defined differently along the corresponding edge. We circumvent these problems by utilizing a doubly-fine grid, i.e., a grid with half the spacing of the original grid. Fig- ure 12b 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 deforma- tion gradient component Fi J and every doubly-fine cell there is precisely one Fi J -value, as specified in Fig. 12a, located on the boundary of the cell. We associate this value to the doubly-fine cell. Thus, a particular value Fi J is distributed to the 4 (in 2D) or 8 (in 3D) adjacent doubly-fine cells, com- pare Fig. 13 for an illustration. The stress components are distributed similarly to the deformation gradients, i.e., the staggering of Fig. 13, is also present for the stresses. With this assignment, to any doubly-fine grid cell all 4 (in 2D) or 9 (in 3D) 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 Fig. 14. We suppose that the constitutive material law is given on the original grid, i.e., each cell is associatedwith onematerial. 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 Pab[ξ ] = 1 8 ∑ l∈L ( Pab l [ξ ] ( Fab l [ξ ] )) ab , (A4) (a) (b) Fig. 13 Placement of deformation gradient variableswithin the doubly- fine material grid: a Location of F11 (shaded in orange). b Location of F12 (shaded in orange). The variables’ grid is depicted by the dashed black lines with ξ = [i, j, k] (A5) L = {l = [l1, l2, l3] with l1, l2, l3 ∈ {0, 1}} (A6) and Paa l [ξ ] = P[ξ ], a ∈ {1, 2, 3}, Pab l [ξ ] = Pba l [ξ ] = P[ξ + l+a,+b], a < b 123 210 Computational Mechanics (2023) 71:191–212 Ta bl e 5 C om po si te bo xe ls fo r m ec ha ni cs — C he at sh ee t Fi ni te st ra in th eo ry In fin ite si m al st ra in th eo ry H ad am ar d ju m p � F � S e = F + − F − = ( 1 c + + 1 c − ) a ⊗ N �ε � S e = ε + − ε − = ( 1 c + + 1 c − ) a ⊗s N D ef or m at io n F ± = F � ± 1 c ± (a ⊗ N ) ε ± = ε � ± 1 c ± (a ⊗s N ) T ra ct io n ba la nc e � P � S e N = 0 �σ � S e N = 0 H es si an fo r N R � f = N ·( A + c + + A − c − ) N � f = N ·( C + c + + C − c − ) N U pd at e of a a[ k] = a[ k− 1] − ( � [k− 1] f ) − 1 f( a[ k− 1] ) a = � −1 f N · (C −− C +) ε � St re ss P � = c + P ++ c − P − σ � = c + σ ++ c − σ − Ta ng en tm od ul us A � = A v − δ A ( N ⊗ � −1 f ⊗ N ) δ A C � = C v − δ C ( N ⊗ � −1 f ⊗ N ) δ C N ot at io n: ± - m at er ia lp ha se s + an d − in co m po si te bo xe l c - vo lu m e fr ac tio n N - no rm al or ie nt at io n of th e la m in at e a - gr ad ie nt ju m p ve ct or a ∈R d F - de fo rm at io n gr ad ie nt ε - in fin ite si m al st ra in P - fir st Pi ol a– K ir ch of f st re ss σ - in fin ite si m al st re ss � f - N ew to n– R ap hs on Ja co bi an � ∈R d ×d A - fo ur th or de r co ns tit ut iv e ta ng en t ∂ P ∂ F = ∂ 2 W ∂ F ⊗ ∂ F C - fo ur th or de r co ns tit ut iv e ta ng en t ∂ σ ∂ ε = ∂ 2 W ∂ ε ⊗ ∂ ε 123 Computational Mechanics (2023) 71:191–212 211 Fig. 14 Placement of the deformation gradients and displacement vari- ables in 3D for l±1,±2 = l±2,±1 = [±l1,±l2, 0] l±1,±3 = l±3,±1 = [±l1, 0,±l3] l±2,±3 = l±3,±2 = [0,±l2,±l3] as well as Faa l [ξ ] = ⎛ ⎜⎝ F11[ξ ] F+1,+2 12,l [ξ ] F+1,+3 13,l [ξ ] F+2,+1 21,l [ξ ] F22[ξ ] F+2,+3 23,l [ξ ] F+3,+1 31,l [ξ ] F+3,+2 32,l [ξ ] F33[ξ ] ⎞ ⎟⎠ , F12 l [ξ ] = ⎛ ⎜⎝ F−1,−2 11,l [ξ ] F12[ξ ] F−2,+3 13,l [ξ ] F21[ξ ] F−1,−2 22,l [ξ ] F−1,+3 23,l [ξ ] F−2,+3 31,l [ξ ] F−1,+3 32,l [ξ ] F−1,−2 33,l [ξ ] ⎞ ⎟⎠ , F21 l [ξ ] = F12 l [ξ ], F13 l [ξ ] = ⎛ ⎜⎝ F−1,−3 11,l [ξ ] F+2,−3 12,l [ξ ] F13[ξ ] F+2,−3 21,l [ξ ] F−1,−3 22,l [ξ ] F−1,+2 23,l [ξ ] F31[ξ ] F−1,+2 32,l [ξ ] F−1,−3 33,l [ξ ] ⎞ ⎟⎠ , F31 l [ξ ] = F13 l [ξ ], F23 l [ξ ] = ⎛ ⎜⎝ F−2,−3 11,l [ξ ] F+1,−3 12,l [ξ ] F+1,−2 13,l [ξ ] F+1,−3 21,l [ξ ] F−2,−3 22,l [ξ ] F23[ξ ] F+1,−2 31,l [ξ ] F32[ξ ] F−2,−3 33,l [ξ ] ⎞ ⎟⎠ , F32 l [ξ ] = F23 l [ξ ], where F±a,±b cd,l [ξ ] = Fcd [ξ + l±a,±b], a �= b. References 1. Arbenz P, van Lenthe GH, Mennel U, Müller R, Sala M (2008) A scalable multi-level preconditioner for matrix-free μ-finite ele- ment analysis of human bone structures. Int J Numer Methods Eng 73(7):927–947. https://doi.org/10.1002/nme.2101 2. Arbenz P, Flaig C, Kellenberger D (2014) Bone structure analysis on multiple GPGPUs. J Parallel Distrib Comput 74:2941–2950. https://doi.org/10.1016/j.jpdc.2014.06.014 3. Kanit T, Forest S, Galliet I, Mounoury V, Jeulin D (2003) Determination of the size of the representative volume element for random composites: statistical and numerical approach. Int J Solids Struct 40(13–14):3647–3679. https://doi.org/10.1016/ S0020-7683(03)00143-4 4. Andrä H, Combaret N, Dvorkin J, Glatt E, Han J, Kabel M, Keehm Y, Krzikalla F, Lee M, Madonna C, Marsh M, Mukerji T, Saenger EH, Sain R, Saxena N, Ricker S, Wiegmann A, Zhan X (2013) Digital rock physics benchmarks - Part I: imaging and segmen- tation. Comput Geosci 50:25–32. https://doi.org/10.1016/j.cageo. 2012.09.005 5. Andrä H, Combaret N, Dvorkin J, Glatt E, Han J, Kabel M, Keehm Y, Krzikalla F, Lee M, Madonna C, Marsh M, Mukerji T, Saenger EH, Sain R, Saxena N, Ricker S, Wiegmann A, Zhan X (2013) Digital rock physics benchmarks - Part II: computing effec- tive properties. Comput Geosci 50:33–43. https://doi.org/10.1016/ j.cageo.2012.09.008 6. Moulinec H, Suquet P (1994) A fast numerical method for comput- ing the linear and nonlinear mechanical properties of composites. Comptes rendus de l’Académie des sciences. Série II, Mécanique, physique, Chimie, astronomie 318(11), 1417–1423 7. Moulinec H, Suquet P (1998) A numerical method for comput- ing the overall response of nonlinear composites with complex microstructure. Comput Methods Appl Mech Eng 157(1–2):69– 94. https://doi.org/10.1016/s0045-7825(97)00218-1 8. Michel JC, Moulinec H, Suquet P (2001) A computational scheme for linear and non-linear composites with arbitrary phase contrast. Int J Numer Meth Eng 52(12):139–160. https://doi.org/10.1002/ nme.275 9. Willot F (2015) Fourier-based schemes for computing themechani- cal response of compositeswith accurate local fields.ComptesRen- dus Mécanique 343(3):232–245. https://doi.org/10.1016/j.crme. 2014.12.005 10. Schneider M, Ospald F, Kabel M (2016) Computational homoge- nization of elasticity on a staggered grid. Int J Numer Meth Eng 105(9):693–720. https://doi.org/10.1002/nme.5008 11. Leuschner M, Fritzen F (2018) Fourier-accelerated nodal solvers (FANS) for homogenization problems. Comput Mech 62(3):359– 392. https://doi.org/10.1007/s00466-017-1501-5 12. Schneider M, Merkert D, Kabel M (2017) Fft-based homogeniza- tion for microstructures discretized by linear hexahedral elements. Int J Numer Meth Eng 109(10):1461–1489. https://doi.org/10. 1002/nme.5336 13. Vondřejc J, Zeman J, Marek I (2014) An FFT-based Galerkin method for homogenization of periodic media. ComputMath Appl 68:156–173. https://doi.org/10.1016/j.camwa.2014.05.014 14. Gélébart Lionel, Ouaki Franck (2015) Filtering material proper- ties to improve FFT-based methods for numerical homogenization. J Comput Phys 294:90–95. https://doi.org/10.1016/j.jcp.2015.03. 048 15. Kabel M,Merkert D, Schneider M (2015) Use of composite voxels in FFT-based homogenization. Comput Methods Appl Mech Eng 294:168–188. https://doi.org/10.1016/j.cma.2015.06.003 16. Kabel M, Ospald F, Schneider M (2016) A model order reduction method for computational homogenization at finite strains on reg- ular grids using hyperelastic laminates to approximate interfaces. 123 https://doi.org/10.1002/nme.2101 https://doi.org/10.1016/j.jpdc.2014.06.014 https://doi.org/10.1016/S0020-7683(03)00143-4 https://doi.org/10.1016/S0020-7683(03)00143-4 https://doi.org/10.1016/j.cageo.2012.09.005 https://doi.org/10.1016/j.cageo.2012.09.005 https://doi.org/10.1016/j.cageo.2012.09.008 https://doi.org/10.1016/j.cageo.2012.09.008 https://doi.org/10.1016/s0045-7825(97)00218-1 https://doi.org/10.1002/nme.275 https://doi.org/10.1002/nme.275 https://doi.org/10.1016/j.crme.2014.12.005 https://doi.org/10.1016/j.crme.2014.12.005 https://doi.org/10.1002/nme.5008 https://doi.org/10.1007/s00466-017-1501-5 https://doi.org/10.1002/nme.5336 https://doi.org/10.1002/nme.5336 https://doi.org/10.1016/j.camwa.2014.05.014 https://doi.org/10.1016/j.jcp.2015.03.048 https://doi.org/10.1016/j.jcp.2015.03.048 https://doi.org/10.1016/j.cma.2015.06.003 212 Computational Mechanics (2023) 71:191–212 Comput Methods Appl Mech Eng 309:476–496. https://doi.org/ 10.1016/j.cma.2016.06.021 17. Kabel M, Fink A, Schneider M (2017) The composite voxel tech- nique for inelastic problems. Comput Methods Appl Mech Eng 322:396–418. https://doi.org/10.1016/j.cma.2017.04.025 18. Uchic MD, Groeber MA, Dimiduk DM, Simmons JP (2006) 3D microstructural characterization of nickel superalloys via serial- sectioning using a dual beam FIB-SEM. Scripta Mater 55(1):23– 28. https://doi.org/10.1016/j.scriptamat.2006.02.039 19. Fliegener S, Luke M, Gumbsch P (2014) 3D microstructure mod- eling of long fiber reinforced thermoplastics. Compos Sci Technol 104:136–145. https://doi.org/10.1016/j.compscitech.2014.09.009 20. Yvonnet J, Bonnet G (2014) A consistent nonlocal scheme based on filters for the homogenization of heterogeneous linear materi- als with non-separated scales. Int J Solids Struct 51(1):196–209. https://doi.org/10.1016/j.ijsolstr.2013.09.023 21. Jänicke R, Diebels S, Sehlhorst H-G, Düster A (2009) Two-scale modelling of micromorphic continua. Continuum Mech Thermo- dyn 21(4):297–315. https://doi.org/10.1007/s00161-009-0114-4 22. Suquet P (1985) Local and global aspects in the mathematical the- ory of plasticity. Plasticity today, 279–309 23. Feyel F (1999) Multiscale FE2 elastoviscoplastic analysis of com- posite structures. Comput Mater Sci 16(1–4):344–354. https://doi. org/10.1016/S0927-0256(99)00077-4 24. Kröner E (1977) Bounds for effective elastic moduli of disordered materials. J Mech Phys Solids 25(2):137–155. https://doi.org/10. 1016/0022-5096(77)90009-6 25. Zeller R, Dederichs PH (1973) Elastic constants of polycrystals. Physica Status Solidi (B) 55(2):831–842. https://doi.org/10.1002/ pssb.2220550241 26. Cooley JW, Tukey JW (1965) An algorithm for the machine calcu- lation of complex fourier series. AMS Math Comput 19(90):297– 301 27. Lahellec N, Michel JC, Moulinec H, Suquet P (2003) Anal- ysis of inhomogeneous materials at large strains using fast fourier transforms 108:247–258. https://doi.org/10.1007/978-94- 017-0297-3_22 28. Eisenlohr P, Diehl M, Lebensohn RA, Roters F (2013) A spectral method solution to crystal elasto-viscoplasticity at finite strains. Int J Plast 46:37–53. https://doi.org/10.1016/j.ijplas.2012.09.012 29. Kabel M, Böhlke T, Schneider M (2014) Efficient fixed point and Newton-Krylov solvers for FFT-based homogenization of elasticity at large deformations. ComputMech 54(6):1497–1514. https://doi. org/10.1007/s00466-014-1071-8 30. Vinogradov V, Milton GW (2008) An accelerated FFT algorithm for thermoelastic and non-linear composites. Int JNumerMethEng 76(11):1678–1695. https://doi.org/10.1002/nme.2375 31. Gélébart L,Mondon-CancelR (2013)Non-linear extension of FFT- based methods accelerated by conjugate gradients to evaluate the mechanical behavior of composite materials. Comput Mater Sci 77:430–439. https://doi.org/10.1016/j.commatsci.2013.04.046 32. Schneider M (2020) A dynamical view of nonlinear conjugate gradient methods with applications to FFT-based computational micromechanics.ComputMech66(1):239–257. https://doi.org/10. 1007/s00466-020-01849-7 33. Schneider M (2021) A review of nonlinear FFT-based computa- tional homogenization methods. Acta Mech 232(6):2051–2100. https://doi.org/10.1007/s00707-021-02962-1 34. Michel JC, Moulinec H, Suquet P (1999) Effective properties of composite materials with periodic microstructure: a computational approach. Comput Methods Appl Mech Eng 172(1–4):109–143. https://doi.org/10.1016/S0045-7825(98)00227-8 35. Leuschner M, Fritzen F (2018) Fourier-accelerated nodal solvers (FANS) for homogenization problems. Comput Mech 62(3):359– 392. https://doi.org/10.1007/s00466-017-1501-5 36. Ospald F, Schneider M, Kabel M (2015) Computational homog- enization of elasticity at large deformations on a staggered grid. In: Conference proceedings of the YIC GACM 2015, pp. 178–191. https://publications.rwth-aachen.de/record/480970 37. Merkert D, Andrä H, Kabel M, Schneider M, Simeon B (2015) An efficient algorithm to include sub-voxel data in FFT-based homog- enization for heat conductivity 105:267–279. https://doi.org/10. 1007/978-3-319-22997-3_16 38. Milton GW (2002) The theory of composites 39. Fritzen F. https://github.com/DataAnalyticsEngineering/ComBo Normal 40. The HDF Group (2022). https://www.hdfgroup.org/ 41. Doll S, Schweizerhof K (2000) On the development of volumetric strain energy functions. J Appl Mech 67(1):17–21. https://doi.org/ 10.1115/1.321146 42. Loehnert S, Mueller-Hoeppe DS, Wriggers P (2011) 3D corrected XFEM approach and extension to finite deformation theory. Int J Numer Meth Eng 86(4–5):431–452. https://doi.org/10.1002/nme. 3045 43. Chen Y, Gélébart L, Marano A, Marrow J (2021) FFT phase-field model combined with cohesive composite voxels for fracture of composite materials with interfaces. Comput Mech 68:433–457. https://doi.org/10.1007/s00466-021-02041-1 Publisher’s Note Springer Nature remains neutral with regard to juris- dictional claims in published maps and institutional affiliations. 123 https://doi.org/10.1016/j.cma.2016.06.021 https://doi.org/10.1016/j.cma.2016.06.021 https://doi.org/10.1016/j.cma.2017.04.025 https://doi.org/10.1016/j.scriptamat.2006.02.039 https://doi.org/10.1016/j.compscitech.2014.09.009 https://doi.org/10.1016/j.ijsolstr.2013.09.023 https://doi.org/10.1007/s00161-009-0114-4 https://doi.org/10.1016/S0927-0256(99)00077-4 https://doi.org/10.1016/S0927-0256(99)00077-4 https://doi.org/10.1016/0022-5096(77)90009-6 https://doi.org/10.1016/0022-5096(77)90009-6 https://doi.org/10.1002/pssb.2220550241 https://doi.org/10.1002/pssb.2220550241 https://doi.org/10.1007/978-94-017-0297-3_22 https://doi.org/10.1007/978-94-017-0297-3_22 https://doi.org/10.1016/j.ijplas.2012.09.012 https://doi.org/10.1007/s00466-014-1071-8 https://doi.org/10.1007/s00466-014-1071-8 https://doi.org/10.1002/nme.2375 https://doi.org/10.1016/j.commatsci.2013.04.046 https://doi.org/10.1007/s00466-020-01849-7 https://doi.org/10.1007/s00466-020-01849-7 https://doi.org/10.1007/s00707-021-02962-1 https://doi.org/10.1016/S0045-7825(98)00227-8 https://doi.org/10.1007/s00466-017-1501-5 https://publications.rwth-aachen.de/record/480970 https://doi.org/10.1007/978-3-319-22997-3_16 https://doi.org/10.1007/978-3-319-22997-3_16 https://github.com/DataAnalyticsEngineering/ComBoNormal https://github.com/DataAnalyticsEngineering/ComBoNormal https://www.hdfgroup.org/ https://doi.org/10.1115/1.321146 https://doi.org/10.1115/1.321146 https://doi.org/10.1002/nme.3045 https://doi.org/10.1002/nme.3045 https://doi.org/10.1007/s00466-021-02041-1 FFT-based homogenization at finite strains using composite boxels (ComBo) Abstract 1 Introduction 1.1 Homogenization in engineering 1.2 FFT-based solvers 1.3 Composite voxel technique 1.4 Outline 1.5 Notation 2 Homogenization problem 3 FFT-based homogenization 4 Composites voxels 4.1 Hadamard jump conditions for finite strain kinematics 4.1.1 Algorithmic implementation 4.2 Robust algorithm via selective back-projection 5 Identification of the normal vector for composite voxels and boxels 5.1 Existing procedures and complications 5.2 Interface indication via a discrete Laplacian 6 Numerical examples 6.1 Material models and loading conditions 6.2 Spherical inclusion 6.3 Composite boxels and local solution field quality 6.4 Selective back-projection in practice 6.5 Short fiber reinforced microstructures 7 Résumé 7.1 Summary 7.2 Discussion Supplementary information Acknowledgements Appendix A A finite difference discretization on a staggered grid References