Mixed strain/stress gradient loadings for FFT-based computational homogenization methods

In this article, the Lippmann–Schwinger equation for nonlinear elasticity at small-strains is extended by mixed strain/stress gradient loadings. Such problems occur frequently, for instance when validating computational results with three-point bending tests, where the strain in the bending direction varies linearly over the thickness of the sample. To control all components of the effective strain/stress gradient the periodic boundary conditions are combined with constraints that enforce the periodically deformed boundary to approximate the kinematically fully prescribed boundary in an average sense. The resulting fixed point and Fletcher–Reeves algorithms preserve the positive characteristics of existing FFT-algorithms, like low memory consumption and extraordinary computational speed. The accuracy and power of the proposed methods is demonstrated with a series of numerical examples, including continuous fiber reinforced laminate materials.


Introduction
The digital design of thin walled lightweight components made of carbon-fiber-reinforced plastics (CFRP) needs accurately calibrated material models on so-called coupon tests [2]. Due to the high in-plane stiffness of the CFRP, tensile measurements turn out to be very difficult, especially the clamping method plays a fundamental role [28]. Therefore, the effective mechanical material parameters are typically measured by three-point or four-point bending tests [3,4,38]. The same holds true for materials like concrete, which have a low tensile strength [1,24].
The measured plate bending stiffness can be predicted for laminate structures made of unidirectional fiber-reinforced laminas, see Fig. 1, very well by using the classical laminate theory [26] in a two-step approach.
1. Calculate the effective lamina stiffness based on the elastic parameters of the fibers and matrix material and the vol-If the lamina exhibits a more complex geometry, e.g., due to fiber waviness, the first step can be replaced by numerical approaches. Especially, the FFT-based homogenization method of Moulinec-Suquet [20,21] has proven to be a powerful tool for the computation of effective mechanical properties of micro-heterogeneous materials. An overview of the improvements of this method is given in [31].
When nonlinear effects enter the stage, a direct simulation of the bending is unavoidable. In the early contribution of Nguyen et al. [22], the authors obtained the force and moment resultants for plates by adding a linear function (with zero mean value) to the ansatz of the strain field and derived a Green's operator for mixed periodic and stress-free boundary conditions. Recently, Gélébart [7] extended the first idea to torsional loadings of beams and applied stress-free boundary conditions by symmetrically extending the plate resp. beam by pore space instead of using the special Green's operator.  In spite of the significant progress accomplished in his contribution, the loading direction still has to be grid-aligned. For practical applications, however, it is often useful to apply a bending in an arbitrary direction without using defective geometrical operations, i.e., rotation of the original unit cell and cutting out a smaller unit cell for computations. From a computational perspective, the extension by pore space typically decreases the convergence speed and enforces the usage of finite element [16,34,42] or finite difference discretizations [33]. Even more severe is the limitation that only the 9 components of the strain gradient which control bending and torsion (see Fig. 3) can be prescribed, and thus direct integration of stress gradient loadings is impossible.
Following Kouznetsova et al. [15] we combine the periodic boundary condition for the displacement fluctuations with kinematic constraints to resolve this issue (see Sect. 2). These constraints enforce the periodically deformed boundary to approximate the kinematically fully prescribed boundary in an average sense and can be automatically fulfilled if the unit cell has certain symmetries with respect to the applied strain gradient loading. E.g., for pure bending loadings it would be sufficient to use PMUBC boundary conditions [10] to obtain zero entries for the 9 additional components of the strain gradient. In the absence of such symmetries, the additional constraints can even strongly influence the solution of pure tensile tests, as we will show with an analytical solution for a two-phase laminate.
Inspired by the higher order homogenization [37,[43][44][45] and multiscale second-order computational homogenization [6,15,40] the kinematic constraints are formulated in terms of strain gradient measures. By transferring the previously developed framework for arbitrary mixed strain/stress boundary conditions [12] to mixed strain/stress gradients in Sect. 3, we derive a Lippmann-Schwinger equation for small strain elasticity with mixed first order boundary conditions.
The algorithms for solving the corresponding fixed point iteration are discussed in Sect. 4. In Sect. 5 we first validate our method with the analytical solution derived in Sect. 2, and then compare it for linear elastic material behavior with the classical laminate theory and numerical results of Nguyen et al. [22] and Gélébart [7] for grid-aligned bending. Afterwards, we extend this comparison to arbitrary loading directions. Finally, we investigate the linear and nonlinear behavior of the 2 mm thick (− 45 • /0 • / + 45 • /0 • ) s 1 laminate shown in Fig. 1.

FFT-based homogenization with strain gradient boundary conditions
Before introducing first order boundary conditions, we will shortly repeat the definition of zero order boundary conditions for FFT-based homogenization [20,21]. For prescribed macroscopic symmetric strain E ∈ Sym 3 and given stress/strain relation σ (ε), where we suppress the x-dependence for notational clarity, we seek a periodic displacement fluctuation u : Y → R 3 on the unit cell Y = [−L 1 /2, L 1 /2] × [−L 2 /2, L 2 /2] × [−L 3 /2, L 3 /2] satisfying the equilibrium equation where ∇ s u = 1 2 (∇u + ∇u T ) denotes the strain due to the displacement fluctuation.
The volumetric components of the macroscopic strain E can be used to apply tensile loads, while the diagonal components allow shear loads to be imposed. For linear elastic material behavior, these load cases are used to compute the effective linear elastic stiffness of the material.

Strain gradient boundary conditions
A naive extension of the equilibrium equation (1) with first order boundary conditions may be obtained by modifying the corresponding optimization problem for the stress potential w, where the subscript # denotes function spaces with vanishing mean value [11]. By adding a linear function E ∇ · x to the argument of the stress potential resp. the definition of ε, the bending response of the material can be computed by solving 1 The lamination scheme of a laminate is denoted following Reddy [26] by (α/β/γ /δ/. . .), where α is the orientation of the first ply, β is the orientation of the second ply, and so on. The plies are counted in the positive thickness direction. This notation also implies that all layers are of the same thickness and made of the same material. For a symmetric laminate, the upper half through the laminate thickness is a mirror image of the lower half. The laminate shown in Fig. 1 is shortly denoted by (− 45 Its solution satisfies the equilibrium equation With E ∇ ∈ Sym 3 3 we want to prescribes the macroscopic gradient of ε. The relevant components of E ∇ for bending loadings are and for torsional loadings These loadings are visualized in Fig. 3. Similar to zero order boundary conditions, where the macroscopic strain E is determined by volume averaging of ε, we want an easy to compute measure for the macroscopic strain gradient E ∇ . Inspired by the work of Kouznetsova et al. [15, equation (A4)] we define the effective bending resp. torsion of the unit cell Y as with component wise division and 1 ∇ ∈ Sym 3 3 denoting the constant strain gradient tensor with only ones, i.e., 1 ∇ i jk = 1 for i, j, k = 1, 2, 3.
By using integration by parts we can show that the following relation holds where · Y + k denotes the average value over Y + k = x ∈ Y : x k = L k /2. Consequently, only the components of the macroscopic strain gradient for bending (5) and for torsion (6) directly coincide with the definition (7) for periodic displacement fields u.

Lippmann-Schwinger reformulation for strain gradient boundary conditions
Definingε = E +∇ s u as well asσ = σ (ε+ E ∇ ·x) allows us to interpret the problem (3) as zeroth-order homogenization where 3×3 denotes the Lippmann-Schwinger operator and G 0 : 3 the solution operator of the linear reference problem, which associates to a right-hand side f the solution of the variational equation 3 . Using the definitions forε andσ we obtain If C 0 is isotropic, i.e., C 0 i jkl = μ 0 (δ ik δ jl + δ il δ jk ) + λ 0 δ ik δ jl with Lamé's moduli λ 0 and μ 0 , this equation can be simplified. Let u 0 3 be the periodic displacement fluctuation defined by and let f 0 and Therefore, a fixed-point iteration based on this extended Lippmann-Schwinger equation (19) allows to prescribe only the components of E ∇ 0 visualized in Fig. 3. This confirms the observation of Gélébart [7], that "...due to the use of periodic boundary conditions, among the 18 strain gradient components, only 9 can be really prescribed.". The remaining 9 components are an outcome of the equilibrium equation (4).
In the following Sect. 2.3 we will heavily rely on Eq. (8) to incorporate additional constraints into our extended Lippmann-Schwinger equation (19) that will enforce the equality of ε ⊗ x ∇ Y and E ∇ for all components.

Kinematically constrained extended Lippmann-Schwinger equation
Without being able to to subject the unit cell to the full gradient E ∇ , it is impossible to incorporate stress gradient loadings into the extended Lippmann-Schwinger equation (19). In other words, more restrictive boundary conditions are needed for mixed strain/stress gradient loadings. Kouznetsova et al. introduced the concept of "generalized" periodicity, which adds the missing 9 constraints [15, equation (25)] Then the deformed boundary approximates the kinematically fully prescribed boundary in an average sense. Therefore, in the following we call the generalized periodic boundary conditions of Kouznetsova et al. more expressively kinematically constrained periodic. According to (8) the kinematic constraints (24) are equivalent to and therefore all components of the macroscopic strain gradient can be obtained by volume averaging of the local strain field The (linear) constrained optimization problem for the bending response under kinematically constrained periodicity condition reads and has the Lagrangian saddle-point formulation By endowing the space H 1 # (Y ) 3 with the Korn-type inner product we can compute the KKT condition as Therefore, the alternating gradient descent ascent method (Alt-GDA) [47] for the primal-dual pair (u, Λ ∇ ) reads for a sequence of positive step-sizes s n and t n . To turn Λ ∇ into a stress gradient we endow additionally Sym 3 3 with the inner product where Y ∇ is the six-order tensor that scales every entry of a third order tensor with the corresponding entry of (1 ∇ · x) ⊗ x Y . Then the gradient of L(u, Λ ∇ ) with respect to Λ ∇ is changed to C 0 : ∇ s u ⊗ x ∇ Y and the Alt-GDA algorithm reads For t n = 1 and s n = 1 we obtain using ∇ s u n = Γ 0 : C 0 : ∇ s u n for the strain fields the following GDA algorithm This can be further simplified to or equivalently This can be transformed into Therefore, the limit strain field ε satisfies the kinematically constrained extended Lippmann-Schwinger equation and its corresponding Lagrange parameters can be computed by

Kinematically constrained periodic boundary conditions
To better understand the effect of the kinematic constraints let us consider the two-phase laminate with thickness h = L 3 shown in Fig. 4 under pure tensile loading in x 1 -direction, i.e.
The two phases of the laminate are isotropic linear elastic with with Lamé's moduli λ ± and μ ± . According to Appendix A of [13] the solution of the unconstrained problem (2), i.e. for periodic boundary conditions, has the phase wise constant strain

Fig. 4 Two-phase laminate
with the rank-one jump given by This solution has the in general non vanishing effective strain gradient component To obtain the solution for kinematically constrained boundary conditions (27) we need to start with a phase wise affine liner ansatz for the strain of the solution. By enforcing the effective strain and strain gradient as well as the continuity of the normal stress at the interface, i.e.
three of the four parameters of the ansatz can be eliminated and we are left with an one dimensional minimization problem for the elastic energy. Simple calculus leads to the following solution ε 11 -Periodic ε 11 -Kinematic. constr. periodic ε 33 -Periodic ε 33 -Kinematic. constr. periodic Due to the corresponding stress does not satisfy the equilibrium equation (1) for γ = 0, i.e. λ − = λ + . In Figs. 5 and 6 the non vanishing components of the strain resp. stress fields are visualized for E + = 1 GPa, E − = 10 GPa, ν ± = 0.3, h = 1mm and E 11 = 10%. Clearly, even at the chosen low phase contrast, the local stress and strain fields differ significantly due to the kinematic constraints. Also the average stresses change notably. Additionally, the solution for kinematically constrained periodic boundary conditions depends on the choice of the center point of the RVE. When choosing a point in the middle of one of the two layers as the center point, the symmetries of the RVE automatically enforce the kinematic constraints and both solutions -with and without kinematic constraints -would coincide.
For the situation at hand, the displacements for kinematically constrained periodic boundary conditions also satisfy zero Dirichlet boundary conditions, compare Fig. 7. Nevertheless, the displacements are completely different inside the RVE, because the kinematically constrained solution does   (1). Put differently, in general it is impossible to satisfy the equilibrium equation (1), periodic boundary conditions and the kinematic constraints all at once. The solution for kinemtically constrained periodic boundary conditions is the energetic optimal deformation that satisfies the periodic boundary conditions and at the same time the kinematic boundary conditions in an average sense.
In Sect. 5 we will use this solution to validate our numerical algorithms and to investigate the local solution quality for different discretizations. Since from an engineering point of view non-symmetric plates that lead to effective bending under tensile loading are not desirable, all other considered examples are symmetric with respect to the center plane. Consequently, the additional constraints would have no effect on the results of (in-plane) tensile simulations and we can use the standard FFT-based algorithms for computing the tensile stiffness then.
Before proceeding with mixed first order boundary conditions, some remarks are in order.

Due to the identity
the constraints (25) are equivalent to which makes the decomposition unique.

The solution of the kinematically constrained extended
Lippmann-Schwinger equation (48) will in general no longer satisfy the equilibrium equation (4). 3. The additional constraints (24) makes it possible to use periodic boundary conditions for the micro scale of a gradient-enhanced FE 2 scheme. Kouznetsova et al. [15] implemented these constraints for the displacements. 4. For displacement based FFT-based schemes [16,18] the constraints can be enforced by using the surface integrals (24) and a kinematically constrained extended Lippmann-Schwinger equation for the displacements similar to (48) can be derived. 5. When the Nyquist frequencies are set to zero [41], the relation (18) is only approximately correct but the resulting linear difference is removed due to the additional constraints (25). 6. The above formulas for E ∇ 0 and E ∇ ⊥ 0 simplify when the reference material is chosen as a scalar multiple of the identity [11, Section 3.1]. Nevertheless, the derivation continues to hold for general isotropic C 0 .

FFT-based homogenization with mixed first order boundary conditions
Thanks to the developments of the previous section, i.e. the kinematically constrained extended Lippmann-Schwinger equation (48), the full gradient E ∇ can be prescribed. We now have the tools in hand to integrate stress gradient boundary conditions, which we will present below, and thus derive a kinematically constrained extended Lippmann-Schwinger equation for mixed first-order boundary conditions.

Stress gradient boundary conditions
Similar to the strain gradient boundary conditions, we must first introduce a measure of the effective stress gradient that is easy to compute.
Using the following extension of the Hill-Mandel energy condition Kouznetsova et al. [15, equation (38)] and Schmidt et al. [29] introduced the effective higher-order stress measure Therefore, we can introduce as an alternative higher order boundary condition the effective stress gradient S ∇ ∈ Sym 3 3 defined by In contrast to the effective strain gradient studied so far, all components of the effective stress gradient can be computed from the local stress field σ for any periodic displacement field u. Before extending the Lippmann-Schwinger equation to mixed first order boundary conditions in the next section, we note that the same extension (68) of the Hill-Mandel lemma is used for higher order homogenization of first strain gradient theories. The macroscopic energy density for elastic material behavior is of the form [37] σ : with L = L 1 = L 2 = L 3 . Consequently, by solving (27) for a suitable sequence of loads E, E ∇ , one can obtain the effective stiffness matrices C 0,0 , C 1,1 and C 0,1 . Keep in mind, that e.g. for the two-phase laminate shown in Fig. 4, the effective stiffness matrix C 0,0 in general does not coincide with the effective stiffness obtained by standard homogenization methods.

Mixed boundary conditions and projectors
To conveniently describe mixed zero order boundary conditions Kabel et al. [12] made use of projectors, i.e., 4-tensors P with minor symmetries which are idempotent P : P : T = P : T for all T ∈ Sym 3 (72) and have major symmetry The first property (72) has the important consequence that every T ∈ Sym 3 can be (uniquely) decomposed in the form Furthermore, T 1 and T 2 can be computed via with the complementary projector Q = I − P, where I is the identity 4−tensor, i.e., I : T = T for all T ∈ R 3×3 . Furthermore, the operators P and Q annihilate each other, i.e., P : Q = 0 and Q : P = 0.
For higher order mixed boundary conditions we additionally need projectors for 3-tensors. For the applications we have in mind it is sufficient to consider only 6-tensors P ∇ that can be represented by three projectors P 1 , P 2 and P 3 for 2-tensors in the following way Using the identity 6−tensor I ∇ , i.e., I ∇ . . . T ∇ = T ∇ for all T ∇ ∈ R 3×3×3 , we obtain the complementary projector Q ∇ = I ∇ − P ∇ . With these preparations in hand, we can neatly formulate mixed boundary conditions for the kinematically constrained extended Lippmann-Schwinger equation (48). Given projectors P, P ∇ and some E, S ∈ Sym 3 as well as E ∇ , S ∇ ∈ Sym 3 3 we seek periodic u : Y → R 3 satisfying the optimization problem (27) and the boundary conditions where ε Y and σ (ε) Y are denoting the average of the strain resp. stress field over Y . Similarly ε ⊗ x ∇ Y is the effective bending resp. torsion and σ (ε) ⊗ x ∇ Y is the effective stress gradient defined in Eq. (7) resp. (70).
Dimension counting shows that (78) and (79) are overdetermined if E and S or E ∇ and S ∇ are unconstrained. Due to the mutual annihilation property (76), a necessary condition for (78) and (79) to be reasonable are the four constraints rendering the optimization problem (27)  1. If P = I and P ∇ = I ∇ , the boundary conditions are equivalent to the constraints ε Y = E and ε ⊗ x ∇ Y = E ∇ . Using for example the loading (in Voigt notation) allows to prescribe the bending visualized in the top left of Fig. 3. 2. By rotating these loadings E and E ∇ by 45 • in the x 1 -x 2 plane, the bending shown in Fig. 8 can be applied. 3. For the projectors P = 0 and P ∇ = 0 the stresses By multiplying the stress gradient with the half of the thickness of a plate resp. the length of a bar, one can easily estimate the maximal occurring tensile stress in the outer layers of a bended plate resp. the maximal shear stress at the end of a twisted bar. 4. Bending of a plate in the x 1 -x 2 plane under plane stress conditions can be performed for P = I−e 3 ⊗e 3 ⊗e 3 ⊗e 3 and P i = I−e i ⊗e i ⊗e i ⊗e i , i = 1, 2, 3, with the loading and In contrast to example 1, the plate can become thinner with this type of boundary conditions due to the Poisson effect. 5. By also rotating the projectors by 45 • in addition to the loads in example 2, the bending shown in Fig. 8 can also be performed under plane stress conditions. A detailed comparison of the results of this simulation setup with laminate theory is presented in Sect. 5.
Similar to the work at hand, Schneider [32] extended the usage of projectors. In his work, he uses projectors to apply mixed zero order boundary conditions in polarization methods.

The kinmatically constrained Lippmann-Schwinger equation for mixed boundary conditions
Given projectors P, P ∇ and prescribed loads E, S ∈ Sym 3 as well as E ∇ , S ∇ ∈ Sym 3 3 satisfying the constraints (80) and (81) the mixed boundary conditions (78) and (79) are incorporated into the optimization problem (27) Following the ideas of Sect. 2.3 we use this constrained optimization problem to derive an extended Lippmann-Schwinger equation that generalizes (48) by mixed boundary conditions. The constrained optimization problem (85) has the Lagrangian saddle-point formulation By endowing the affine linear subspaces P :ε = E of Sym 3 and P ∇ . . .ε ∇ = E ∇ of Sym 3 3 with the inner products and keeping the inner product for the space Sym 3 3 of Lagrange parameters Using the Moore-Penrose pseudoinverse [19,25] M = Q : C 0 : Q † as well as M ∇ = Q ∇ . . . C 0 : Q ∇ † whose components are given by we can compute the KKT condition as Therefore, the alternating gradient descent ascent (Alt-GDA) [47] reads for a sequence of positive step-sizes s n and t n . For t n = 1 and s n = 1 we obtain using ∇ s u n = Γ 0 : C 0 : ∇ s u n for the strain fields the following GDA algorithm Due tō the update of the average strain and its gradient can be written asε The two other equations can be simplified to or equivalently This can be transformed into so that the complete algorithm using the variable τ for the polarization reads Therefore, the limit strain field ε satisfies the kinematicall constrained extended Lippmann-Schwinger equation with mixed first order boundary conditions where Γ 0 : ·⊗x ∇ Y ,⊥ 0 and M ∇ . . . Q ∇ . . . ·⊗x ∇ Y are shorthand notations for the operators and The corresponding Lagrange parameters of the strain field ε can be computed by The equation (122) can be simplified if the projectors P, P ∇ and the reference tensor C 0 commute, i.e., C 0 : P : T = P : C 0 : T , holds for all T ∈ R 3×3 and T ∇ ∈ R 3×3×3 . This is automatically satisfied if C 0 is chosen as a multiple of the identity tensor I. Conversely, (127) holds for all projectors P, P ∇ if and only if C 0 is a multiple of the identity. Then, equation (122) simplifies to where D 0 denotes the inverse of C 0 . An alternative derivation more in the line of [12] can be found in the Appendix 1

FFT-based algorithms
In this section we will derive FFT-based algorithms to solve the extended Lippmann-Schwinger equation (122). The projected gradient descent method [23,27] will use Λ ∇,n+1 ⊥ 0 instead of the last iterate Λ ∇,n ⊥ 0 for the update of ε n in equation (121) and can therefore be seamlessly integrated into existing (nonlinear) conjugate gradient (CG) methods [30,46].

Projected gradient descent
After the derivations of Sect. 3.3 it seems natural to use the gradient descent ascent method [47] for the solution of the Fig. 9 Projected gradient descent method extended Lippmann-Schwinger (122). It is known, however, that even the convergence speed of gradient descent methods is limited [11] and that one should use (nonlinear) CG methods whenever possible [30].
Since it seems to be difficult to integrate the gradient descent ascent algorithm (116)-(121) into a CG method, we decided to use the projected gradient descent (PGD) method [23,27]. This method projects the iterates of the standard gradient descent method for the unconstrained problem onto the (affine linear) subspace of admissible functions satisfying the constraints. From an algorithmic point of view this can be achieved by only changing the update of the strain field in the gradient descent ascent algorithm. More precisely, using Λ ∇,n+1 ⊥ 0 instead of the last iterate Λ ∇,n ⊥ 0 for the update of ε n allows to eliminate the Lagrange parameters While the first part (129)-(132) performs the classical basic scheme of Moulinec-Suquet on the unconstrained extended Lippmann-Schwinger equation, the last step (133) projects the next iterate on an admissible state. Unfortunately, as the visualization in Fig. 9 shows, the projection is in general not orthogonal to the subspace of admissible solutions.  Inverse Fast Fourier transform 14: Strain gradient with constraints (131),(135) Adding linear strain 16: return ε Algorithm 2 CG method for mixed boundary conditions (multiple load steps) g old ← g 14: g ← G 2 15: D ← g gold D − G 16: until Convergence 17: end for 18: return ε Algorithm 3 Fletcher-Reeves method for mixed boundary conditions (multiple load steps) 1: ε ← 0 Initialization of ε 2: for k ∈ {0, . . . , k max } do 3: γ ← 1 4: until Convergence 13: end for 14: return ε The projection can only be interpreted as ascent step in orthogonal direction to the (affine linear) subspace of admissible solutions.
Using the equivalent steps instead of (132) and (133) this method can be implemented as shown in Algorithm 1.
Notice that the operators M = Q : C 0 : Q † and M ∇ = Q ∇ . . . C 0 : Q ∇ † have to be computed only once in preprocessing of each loading step, when C 0 is adapted by setting it to the average of the maximal and minimal (positive) eigenvalue of the tangential stiffness dσ/dε [11]. Similar as described by Kabel et al. [12] for unidirectional tensile tests, the following affine linear extrapolation improves the convergence speed for strain gradient simulations without superimposed stretch

Conjugate gradient methods
In the case of linear material behavior, i.e. σ (ε) = C : ε, the extended Lippmann-Schwinger equation with mixed first order boundary conditions (122) is a linear relation and therefore the conjugate gradient (CG) method can be used. This method was first applied by Zeman et al. [46] for FFT-based homogenization. It needs four instead of one solution vector but can be expected to have a significantly improved convergence behavior, see the detailed comparison of Schneider [30] for standard boundary conditions. The implementation shown in Algorithm 2 is based on the Moulinec-Suquet iterate for computing the matrix vector product on the left hand side of (122), which was previously used as core function of the PGD implementation shown in Algorithm 1.
In the case of nonlinear material behavior we suggest using the Fletcher-Reeves (FR) nonlinear CG [30] for solving (122), see Algorithm 3. The implementation proposed by Schneider [30] can be stabilized by limiting the factor γ γ old , compare line 10 of Algorithm 3. It needs one solution vector less but typically also converges slower than the standard CG method in the linear elastic case. To reduce memory requirements, both the CG and FR methods could also be implemented by using additional displacement vectors instead of additional strain vectors [9,11].
In Sect. 5.4 we will compare the convergence speed of all three methods.

Examples
In the following we will discuss multiple examples of laminate structures. At first, we will use in Sect. 5.1 our analytical solution for the two-phase laminate shown in Fig. 4 under pure tension with the additional constraints to study the influence of the discretization on the local solution quality. After discussing the relation of the force and moment resultant with the effective stress and stress gradient in Sect. 5.2 we will then investigate in Sect. 5.3 the examples of Nguyen et al. [22] and Gélébart [7] for which only grid-aligned loadings are applied. Afterward, we will modify one of their examples and validate our extended Lippmann-Schwinger equation in the linear elastic regime by comparison with the classical laminate theory [26] in Sect. 5.4. Finally, in Sect. 5.5 we will study the effective linear and nonlinear behavior of a multilayer laminate subjected to tensile and bending loadings. All simulations results were obtained with FeelMath [5] which is distributed as part of GeoDict 2 .

Strain gradient boundary conditions
In this section, the analytical solution shown in Figs. 5, 6 and 7 for a two-phase laminate under tensile loading and with suppressed strain gradient is used to validate our method and to investigate the local solution quality.
In the middle column of Fig. 10 we compare the local strain fields for three different resolution and two discretizations. Obviously, the staggered grid discretization [34] exhibits a significantly higher local accuracy than the Hex8R discretization [42]. This is at least partly related to the fact, mentioned earlier, that for the Hex8R discretization the Nyquist frequencies are set to zero [41]. Somewhat surprisingly, the Hex8R discretization does not show comparable solution quality even at doubled resolution. Only the (phasewise) averages are of similar accuracy, compare Table 1.
Since the displacement field is obtained by integrating the strain field according to the formulas derived in [10], we make the same observation for the local displacement field shown in the left column of Fig. 10.
The non-trivial stress components σ 11 and σ 22 are not shown because their errors are only a constant rescaling of the error for σ 33 that is visualized in the right column of Fig. 10. This is due to the fact that the numerical error for ε 11 is equal zero. Therefore, the phase-wise scaling factor of the error coincides with λ ± /(λ ± + 2μ ± ), which for both phases is equal to 3/7.
Despite these observations, we will only use the Hex8R discretization in the following. Otherwise, our results would not be comparable to the benchmark results of Gélébart [7]

Force and moment resultant
To be able to compare our simulation results for a plate thickness h in the x-y-plane with the benchmark results of Nguyen et al. [22] and Gélébart [7] in Sect. 5.3 as well as classical results of the laminate theory [26] in Sects. 5.4 and 5.5 we first discuss the relation between the force resultant N and the moment resultant M with the effective stress σ Y and the effective stress gradient σ ⊗ x ∇ Y of the plate.
Since the components of the force and moment resultants are defined by for i, j ∈ {1, 2} we obtain where To simulate plates with plane stress assumption and to directly obtain coefficients of the extensional stiffness matrix A and bending stiffness matrix B we apply the following loadings in Voigt notation 1. Tension: P = e 1 ⊗ e 1 ⊗ e 1 ⊗ e 1 + e 2 ⊗ e 2 ⊗ e 2 ⊗ e 2 + e 1 ⊗ e 2 ⊗ e 1 ⊗ e 2 and P ∇ 1 = P ∇ 2 = P ∇ 3 = I with and Then -depending on the applied loading E 11 resp. E ∇ 113it is possible to directly compute entries of the extensional resp. bending stiffness matrix from the local stress field 1. Tension: 2. Bending: If Y is symmetrically extended by pore space in thickness direction, as proposed by Gélébart [7], the same formulas apply when h is replaced by the thickness of the extended plate.
The closed-form solution for the effective stiffnesses of the Kirchhoff plate with n layers are given by the following relations [26] where C k denotes the two-dimensional stiffness matrix under plane stress assumption for the kth layer of the plate located between the points x 3 = h k−1 and x 3 = h k in thickness direction. Under the assumption of isotropic material behavior described by Young's modulus E and Poison's ratio ν the stiffness matrix is given by and, more generally, for transversely isotropic material behavior with Young's moduli E 1 , E 2 , in-plane shear modulus G 12 and in-plane Poisson's ratio ν 12 (ν 21 = ν 12 E 2 /E 1 ) The four independent parameters are related to the five independent parameters of a transversely isotropic material in the following way The plate stiffness matrices A, B and D of a rotated plate can be obtained by rotating the stiffness matrices C k for each layer [8].

Grid-aligned loading directions
In the following we compare the results of our approach for grid-aligned loading directions with the benchmark results of Nguyen et al. [22] and Gélébart [7]. The first benchmark example is a laminate plate with moderate phase contrast, see Fig. 11 and Table 2. In this case, the results with and without extension by pore space coincide if the correct thickness is used for the formulas (146) and (147). Therefore, in Tables 3  and 4 we only show the values without extension. Since we do not increase the phase contrast by an extension, our approach converges for the conjugate gradient method shown in Algorithm 2 with fewer iterations and also gives more precise results at the same time. For the second benchmark example of a plate with periodic inclusions and three different sets of material parameters, see Fig. 11 and Table 2, it is not possible to reproduce the FEM reference solution of Nguyen et al. [22] without using an extension by pore space. Mixed boundary conditions [12] only allow to prescribe zero average stress but not zero surface stress as prescribed in Nguyen et al. [22]. The resulting stress fields for the third parameter set are visualized in Fig. 12 for tensile and in Fig. 13 for bending loadings. Consequently, we show in Tables 5, 6, 7, 8, 9 and 10 the results of our approach for the microstrucuture with pore space extension to obtain comparable results under tensile loadings. Following the observation of Gélébart, we used mixed boundary conditions, i.e., P = e 1 ⊗ e 1 ⊗ e 1 ⊗ e 1 + e 2 ⊗ e 2 ⊗ e 2 ⊗ e 2 + e 1 ⊗ e 2 ⊗ e 1 ⊗ e 2 , on the extended geometry to reduce the number of iterations. Nevertheless, for this infinite contrast problem only the conjugate gradient method 2 and the Fletcher-Reeves algorithm 3 converge at a moderate number of iterations. Fig. 11 Microstructures used by Nguyen et al. [22] and Gélébart [7]. Top: Laminate plate. Bottom: Plate with periodic inclusion. Middle: Central part of the plate with periodic inclusions which is first homogenized in the two-step homogenization approach in Sect. 5.4    For these results, we did not use the energy based convergence criterion proposed by Nguyen et al. [22] because it leads to inaccurate solutions when combined to our approach. With our strain based convergence criterion [11] ε new 2 − ε old 2 we obtain in general more accurate solutions than Nguyen et al. [22] and Gélébart [7] especially for the coarsest resolution. On the other hand, we also need significantly more iterations for convergence. Without the extension by pore space, this would not be the case. We will study the convergence behavior and the influence of the extension on the convergence speed in detail for arbitrary loading directions in the next section.

Arbitrary loading direction
To validate the derived Lippmann-Schwinger equation for mixed boundary conditions, we compare in the following the results of our approach with the 2-step homogenization approach described in the introduction of this paper. Contrary to the previous section, we now use mixed (strain/stress) boundary conditions only on the geometries without extension.
As test geometry we use the plate with inclusion from the previous section with the third parameter set, i.e., E inclusion = 10 GPa and E matrix = 1 GPa. Numerically homogenizing the central part of the plate with inclusion shown in the middle of Fig. 11 we obtain the anisotropic stiffness parameters shown in Table 11. In Fig. 14 we compare five different results.
1. The results of the laminate theory using the parameters of Table 11. 2. Numerical results for the laminate structure without extension by pore space. 3. Numerical results for the laminate structure with extension by pore space. 4. Direct numerical results on the plate with inclusion without extension by pore space. 5. Direct numerical results on the plate with inclusion with extension by pore space.
The left part of Fig. 15 shows that the first three results coincide independently of the loading direction. The last two results on the resolved geometry with and without extension coincide only for the bending stiffness but not for the extensional stiffness due to the boundary effects discussed in the previous section. Then the plane strain and zero surface stress boundary conditions lead to a softer response.
The relative difference of the results on the resolved geometry and the two-step approach is below 4% for the extensional and below 8% for the bending stiffness, see right part of Fig. 15. For tensile and bending loadings in fiber direction (y-direction) the five results are even in perfect agreement.
The number of iterations necessary for convergence, shown on the right side of Fig. 14 is robust with respect to the loading direction and the resolution of the geometry, but extending the geometry by pore space significantly decreases the convergence speed.

Multilayer laminate
For representative volume elements the boundary effects, enforcing us to use an extension by pore space in Sect. 5.3, are by definition negligible. Therefore, for our representative microstructure of a multilayer laminate, see Fig. 1, we will only discuss the results obtained without extension. The (-45 • /0 • /+45 • /0 • ) s laminate has a total thickness of 2mm. Each unidirectional layer has a carbon fiber volume content of 40%. The transversely isotropic elastic material parameters of the carbon fibers are shown in Table 12. Since for bending always a combination of tensile and compressive loading occurs, the tension-compression asymmetry of epoxy has to be taken into account by the material model. Therefore, we model the epoxy matrix material by an elasto-plastic material model of Stier et al. [36] that incorporates different yield strengths in tension and compression. The material parameters of Table 13 were calibrated by Ullah et al. [39] according to the measurements of Stier et al. [36] on pure epoxy samples, see Fig. 16.
At first we compare -similar to the previous section -in the first row of Fig. 17 the extensional and bending stiffness predicted by the following three methods. The iteration numbers shown in the second row of Fig. 17 reveal that also for this complicated geometry, the convergence speed is not influenced by the loading direction. As to be expected for the homogenized unidirectional layers the FFT-based homogenization methods needs fewer iterations. The last row of Fig. 17 shows the difference between the predictions of the three approaches. The two-step numerical homogenization approach is closer to the laminate theory, because it uses the same intermediate results. For any loading direction, the relative difference of the extensional and bending stiffness is below 0.15% resp. 0.5%.      Table 10 Bending stiffness obtained for the inclusion case (with extension by pore space) for E inclusion = 10 GPa and E matrix = 1 GPa

Resolution
Analytical Nguyen et al [22] Iter . Gélébart [7] Iter.  Table 11 Anisotropic stiffness of the central part of the plate with inclusions shown in the middle of Fig. 11 for the third parameter set of Table 2 Parameter At this point one could argue, that deriving an extended Lippmann-Schwinger equation is not really helpful, because it only reproduces the results of classical laminate theory and the original version of the Lippmann-Schwinger equation [20,21] would be sufficient to derive the stiffness of the unidirectional layers. When nonlinear effects enter the stage, this is no longer true. Only the direct simulation approach can be applied to physically nonlinear material behavior without any changes. On the left of Fig. 18 16 Calibrated material model for epoxy [39] the stacking order of the unidirectional layer does not influence the force resultant. Therefore, the force resultant for α • and 180 • − α • coincide. For the moment resultant under bending loadings this is not true because the outer layers of the (-45 • /0 • /+45 • /0 • ) s laminate, e.g., the −45 • layers, are submitted to higher loadings. By post-processing of these force and moment resultants and the corresponding evolution of the average plastification, the yield surfaces under tensile and bending loadings shown in Fig. 20 were derived. For tensile loadings the yield stress is the highest in the main fiber direction of 0 • . Under bending loadings the yield surface is reoriented into the direction of the outermost fiber orientation.
The computational effort for tensile and bending simulations could be kept limited thanks to the application of the Fletcher-Reeves algorithm 3. The iteration numbers for convergence of each loading step on the right of Fig. 18 were obtained with affine linear extrapolation (136) between the loading steps. Only for the first iterate, these numbers are independent of the loading direction. Then the plastic evolution changes the convergence behavior significantly. Since bending loadings lead to higher local plastifications than comparable tensile loadings, the iteration numbers are approximately 50% higher for this type of loading.

Conclusion
In this work, we first extended the Lippmann-Schwinger equation to strain gradient boundary conditions. To be able to prescribe the full strain gradient of the solution, it was necessary to introduce additional kinematic constraints. This allowed us to also integrate stress gradient boundary condi-tions and to enforce mixed first order boundary conditions with respect to an arbitrary coordinate system. Thanks to the use of a projected gradient descent method, the kinematic constraints can be easily implemented in existing FFT-based algorithms. This extension makes it possible, for example, to simulate plate bending. In contrast to the simulation of tensile experiments, a volume element can only be representative for  When applied to the test cases with grid-aligned loadings of Nguyen et al. [22], our method gave more accurate predictions at coarse resolutions. By extending the example to anisotropic material behavior, we could validate our approach for arbitrary loading directions with the classical laminate theory. It turned out that the approach of Gélébart [7] to extend the domain by pore space significantly decreases the convergence speed. The difference in the boundary conditions due to the pore space extension, e.g., zero surface stress instead of periodic stress, did not influence the predicted effective response.
Regarding the multilayer laminate example, in the linear elastic regime the predictions of our approach coincide with the classical laminate theory, if the boundary conditions are set properly. By applying the Fletcher-Reeves method [30] we were able to predict the force and moment resultant under arbitrary loading directions at moderate computational effort. Thereby, we could also determine the yield surface of the laminate under tensile and bending loadings.
It should be remarked, that our derivation directly carries over to problems of finite strains. Therefore, our FFT-based solver for the kinematically constrained extended Lippmann-Schwinger equation could replace-similar to previously developed FE-FFT methods [14,35]-the microscale solver in higher order FE 2 methods [15,29,40] at small and finite strains.  Invoking the splitting