Imposing different boundary conditions for thermal computational homogenization problems with FFT‐ and tensor‐train‐based Green's operator methods

To compute the effective properties of random heterogeneous materials, a number of different boundary conditions are used to define the apparent properties on cells of finite size. Typically, depending on the specific boundary condition, different numerical methods are used. The article at hand provides a unified framework for Lippmann–Schwinger solvers in thermal conductivity and Dirichlet (prescribed temperature), Neumann (prescribed normal heat flux) as well as periodic boundary conditions. We focus on the explicit jump finite‐difference discretization and discuss different techniques for computing the discrete Green's operator. These include Fourier‐type methods, that is, based on discrete sine, cosine and Fourier transformations, as well as low‐rank tensor methods, that is, the tensor‐train approximation, whose computational prowess was demonstrated recently. We use the developed computational technology to investigate a number of interesting random materials and assess different microstructure‐generation techniques in terms of their compatibility with the prescribed boundary conditions.


INTRODUCTION
Many materials used in industry show a random heterogeneous microstructure. 1 Computational homogenization approaches permit to simulate components made from materials with a heterogeneous microstructure whose microstructure features are orders of magnitude smaller than the part under consideration by identifying an equivalent effective material model (such as an effective thermal conductivity or an effective stiffness) via a mathematical limiting procedure.Periodic homogenization is applicable to periodic microstructures, whereas stochastic homogenization is required to deal with materials that show a random microstructure.
The effective properties of a random material are computed on a representative volume element, that is, a cell that is statistically typical 2,3 of the whole part and the influence of the boundary conditions on the effective properties is negligible. 2 For periodic homogenization, it is natural to apply periodic boundary conditions for computing the effective properties on the periodic unit cell.The choice of boundary conditions for computing the apparent properties in stochastic homogenization is not that straightforward as any subvolume of a random heterogeneous microstructure is, in general, nonperiodic.
A particularly performing class of numerical methods on regular grids was Moulinec and Suquet, 4,5 which are based on the fast Fourier transform (FFT).These methods are oftentimes less resource-intensive than classic finite element (FE) methods for a number of reasons, for example, by discretizing the domain on a fixed grid and by employing efficient preexisting implementations of the FFT.For FE, applying the field values (Dirichlet boundary condition) or the gradient of the field value (Neumann boundary conditions) on the boundary is natural from the method's formulation whereas periodic boundary conditions require some additional effort. 6The differences between Dirichlet, Neumann, and periodic boundary conditions were studied for FE discretizations. 7FT-based solvers offer significant performance advantages for the computation of apparent properties over FE approaches, but inherently impose periodic boundary conditions due to the periodic nature of the Fourier transform.Under the constraint of periodic boundary conditions, there are two ways to handle random microstructures.The first option is to generate a synthetic microstructure that resembles the nonperiodic microstructure in all relevant measures but is periodic, thereby obtaining a microstructure that aligns with the imposed periodic boundary conditions.The second option is to use the nonperiodic microstructure, called snapshot by Schneider et al. 8,9 and to impose periodic boundary conditions on this part, accepting that the microstructure does not align with the boundary conditions.In this case the statistical homogeneity is disturbed. 8Nonetheless, in case of the "snapshot" microstructures, the use of FFT-based methods with periodic boundary conditions provides a middle-ground between the overestimation of the apparent stiffness by Dirichlet boundary conditions and the underestimation of the stiffness by Neumann boundary conditions, an insight proposed by Hill 2 in the 1960s and corroborated later on. 10,11everal approaches for homogenization using FFT-based methods using nonperiodic boundary conditions requiring additional effort were proposed. 12Neumann boundary conditions were proposed for FFT-based methods by doubling the domain 9 in analogy to existing approaches for solving the Poisson equation. 13Dirichlet boundary conditions were applied using FFT-based methods by using an enlarged domain with an elastic buffer zone, effectively working on a subdomain of an FFT-based solution field. 14Beyond the scope of FE and FFT-based methods Neumann and Dirichlet boundary conditions were applied using the conventional FE discretizations.Khoromskaia et al. used a tensor-train-based 15 approach to apply Neumann and Dirichlet boundary conditions. 16The solution of Poisson's equation, which is algorithmically similar to the cell problem for thermal conductivity, was performed with FFT-based approaches with Dirichlet, Neumann, and periodic boundary conditions. 17,18A possible solution to the similar problem in linear elasticity was introduced by means of the mixed uniform boundary conditions, 19,20 allowing for a handling of nonperiodic boundary conditions, although not in the same manner in which periodic boundary conditions are treated.For wave propagation, Dirichlet boundary conditions were handled using a Schur-complement preconditioner for solving a KKT system. 21or computational homogenization, some caution is advised when using the term "boundary condition."In fact, the type of imposed macroscopic field, that is, the strain or the stress or a mixture of both, is sometimes called the "boundary condition."However, it is well possible to use periodic boundary conditions for the displacement-fluctuation field and impose the average stress, for instance.
A variety of discretization schemes were proposed for FFT-based computational homogenization methods and periodic boundary conditions.In fact, the discretization by Moulinec and Suquet, 4,5 the finite difference 22 and rotated staggered grid [23][24][25] discretization schemes proved to be adequate compromises between numerical stability and the required computational effort. 26Another discretization with a strong similarity to the finite difference discretization suitable to heat transfer problems is the explicit jump (EJ) discretization. 27he pioneering solution method for FFT-based computational homogenization was the so-called basic scheme by Moulinec and Suquet 4 which can be interpreted as a preconditioned gradient descent scheme. 28Based on this insight, improved solution methods were proposed to solve the Lippmann-Schwinger equations, such as the gradient method with adaptive step sizes proposed by Barzilai and Borwein. 29,30Additionally, other solvers, such as the linear 31,32 and nonlinear conjugate gradient method 33 or Newton's method 28,34 were introduced.We refer to the dedicated reviews of FFT-based homogenization [35][36][37][38][39] for further insight into discretization schemes and solution methods.Most of these solvers are formulated with the gradient of the quantity of interest as the solution variable (i.e., the strain or stress field for a mechanical problem or the temperature gradient or heat flux for a thermal problem).In some cases, the handling of boundary conditions can be simplified by choosing solution schemes formulated directly on the quantity of interest, for instance the temperature or displacement field. 26,40,41 key ingredient for the effectiveness of FFT-based methods is the efficient computation of the preconditioner for the gradient descent scheme. 28In case of displacement or temperature based solvers, the preconditioner applies the inverse of a multidimensional Laplacian in Fourier space.The inverse Laplacian is also referred to as Green's operator.This Green's operator may also be applied in real space using a convolution approach 42 or a low-rank method. 16For certain applications, for example, homogeneous finite difference operators, many entries of the array are redundant. 43ow-rank methods work by decomposing a given array into a product of lower-dimensional arrays.For the canonical decomposition, 44,45 the array is decomposed into a sum of outer products of vectors.Although algorithms to approximate arrays by canonical decompositions [46][47][48][49] exist, these algorithms were proven 50 to be intrinsically unstable.Another available decomposition is the Tucker decomposition, 51 which, although it can reduce the memory required, is not free from the curse of dimensionality. 52Any low-dimensional array-like structure (like the finite-difference Laplacian) can be reshaped as a high-dimensional array with very low number of entries per dimension in case the number of grid points per dimension has a prime decomposition with small prime numbers.By selecting the number of entries per dimension of this reordered matrix as two and restricting to original arrays with sizes of a power of two, we employ a low-rank method for high-dimensional arrays, called the tensor train (TT) format 15 also known as the matrix product state in quantum physics 53,54 on low-dimensional arrays occurring in solution schemes for computational homogenization.Reshaping a low-dimensional array into a high-dimensional one and applying the TT decomposition is called the quantics tensor train (QTT) format. 55,56n fact, the finite-difference Laplacian allows for an exact low-rank representation in the QTT format without any loss of accuracy.This approach therefore seems suitable as a method for inverting a high-dimensional Laplacian and to approximate Green's operator for the solution of homogenization problems.This low-rank approximation of Green's operator was studied for stochastic homogenization of small-scale problems. 16,57Such a strategy permits to impose different types of boundary conditions, as proper low-rank representations (in the QTT format) exist for the Laplace operator with the respective boundary conditions.An alternative strategy based on storing low-rank approximations of the involved field was studied by Vondřejc et al. 58 Applying Green's operator in Fourier space brings some advantages that cannot be mitigated by real-space methods.Therefore, it seems desirable to study approaches based on Fourier-related transforms.For instance, Dirichlet and Neumann boundary conditions may be enforced via the discrete sine transform (DST) and the discrete cosine transform (DCT).This approach is standard for solving Poisson's equation, 59,60 but appears to be less used for problems with heterogeneous coefficients.
In this work, we present a unified approach to applying Neumann and Dirichlet boundary conditions when computing the effective conductivity of heterogeneous microstructures.The remainder of this article is structured as follows.In Section 2, we derive a unified formulation of cell problems for all three types of boundary conditions.Using this unified formulation, in Section 3, we present a solution scheme covering all types of boundary conditions.Therefore, we formulate the solution scheme to the cell problem for all boundary conditions by introduction of Green's operator as a preconditioner and derive the corresponding convergence estimates.We recapitulate the accelerated gradient method based on the work of Barzilai and Borwein for temperature-based homogenization.We revisit the EJ discretization and show how boundary conditions are treated in the discrete case using the unified framework from Section 2. Next, we show how to evaluate the corresponding Green's operator for Dirichlet, Neumann, and periodic boundary conditions using two different approaches.The first and more standard approach uses the DST/DCT to replace the FFT for the application of Green's operator with Neumann and Dirichlet boundary conditions, which we will refer to as the FFT-based approach.In addition, we show a TT-based approach to Green's operator which works by exploiting the low-rank TT format to perform a matrix inversion of large matrices via exponential sums.Lastly, in Section 4, we compare both approaches of computing Green's operator in terms of the speed of convergence and computational effort.Then, using the FFT-based approach, we evaluate numerical examples to assess the influence of boundary conditions on computing the effective conductivity of large-scale microstructures.

The effective conductivity
We are concerned with a thermally conducting random material, given in terms of a random field of conductivity tensors in d = 2, 3 dimensions.More precisely, we assume that a statistical ensemble of conductivity fields is given, and we denote by ⟨⋅⟩ the associated ensemble average.We suppose that there are universal positive constants  ± , s.t., with probability one, the random conductivity fields (1) are symmetric and satisfy the inequalities for almost every x ∈ R d and all  ∈ R d .Following the seminal papers of Hill, 2 there is a number of conditions which ensure the emergence of a mathematically well-defined macroscopic behavior of such a random material.For a start, the random field is assumed to be stationary, that is, for any integrable random variable , shifting by a vector z ∈ R d preserves the expectation of the random variable Moreover, we assume the ensemble to be ergodic, that is, for any integrable random variable , we may express the ensemble average by a spatial average over the whole space for almost every realization .Stationarity (4) encodes a form of stochastic translation invariance whereas ergodicity (5) ensures that a single realization provides "rich enough" values to recover an independent sampling.In stochastic homogenization, 61 which presupposes scale separation, 2 the effective thermal conductivity  eff ∈ R d ⊗ R d is a deterministic tensor which describes the average behavior of the random thermally conducting medium responding to the field (1) on large scales.
To describe a conductivity tensor, it is sufficient to prescribe its action on arbitrary negative* temperature gradients  ∈ R d .Thus, we fix such a negative temperature gradient  ∈ R d .Then, the effective conductivity in direction  is defined as where  ∶ R d → R is the (temperature) corrector and the evaluation takes place at an arbitrary location x ∈ R d .The corrector  is a random field which grows at most linearly and solves the static heat equation (in a weak sense) on the whole space R d .It can be shown 62,63 that the corrector exists and is unique under a suitable normalization condition, for example, having vanishing mean on the unit ball.Furthermore, it turns out that the gradient of the corrector ∇ defines a stationary random field.As a consequence of this stationarity, the definition (9) of the effective conductivity does not depend on the point of evaluation x ∈ R d .With this definition (9) at hand, the corrector equation ( 7) permits to deduce that the effective conductivity tensor  eff is symmetric and obeys the bounds (2) The definition ( 9) is rather convenient to deduce these properties, but cumbersome to evaluate.In fact, computing such an ensemble average may come with significant effort.With the help of the ergodicity property (5), it is possible to express the effective conductivity as the volume average over the corrector (7) corresponding to a single realization  (with probability one).Although the formula (9) for the effective conductivity involves only a single realization, it still requires to compute an integral over the entire space R d .Moreover, the integrand entails the corrector  which solves the corrector equation ( 7) on the whole space, as well.

2.2
The apparent conductivity with different boundary conditions

Generalities
To compute the effective conductivity (9) in practice, it is imperative to devise approximations that involve operations on finite domains, typically cubes ] d (10)   with edge length L > 0. We consider a stationary and ergodic proxy random field of conductivities  L .Upon first reading, it is instructive to assume that the original random field  and the proxy field  L are identical.However, in a number of practical situations the extra freedom of working with the field  L is actually critical.For instance, considering periodized ensembles, that is, conductivity fields  L which are periodic on the cell Y L , often lead to higher accuracy than considering cut ensembles, 64 that is, those with  ≡  L .
When working on the cell (10), it is customary to seek appropriate replacements for the full-space corrector (7).However, to give rise to well-posed problems, elliptic partial differential equations need to be supplemented with appropriate boundary conditions.

Dirichlet boundary conditions
For instance, we may define the Dirichlet corrector  D,L ∶ Y L → R as the solution of the boundary-value problem It is not difficult to see that the problem (11) has a unique (weak) solution in a suitable Sobolev space.Then, the apparent Dirichlet conductivity is defined via Notice that the apparent conductivity  app D,L is still a random variable, in contrast to the effective conductivity (9).However, the apparent conductivity (12) shares some salient properties of the effective conductivity.The defining Equation (11)  permits to show that the apparent conductivity (12) is almost surely symmetric and satisfies the bounds Last but not least, under the condition that the proxy field  L converges, in a suitable sense, to the original coefficient field (1), the apparent conductivity (12) converges to the effective conductivity lim

Periodic boundary conditions
Choosing Dirichlet, that is, temperature, boundary conditions is a popular choice for the cell problem (11).However, other boundary conditions are frequently used, as well.For instance, one may consider the periodic corrector  P,L ∶ Y L → R, defined as the solution to the boundary-value problem which encodes the continuity of the temperatures and of the normal fluxes across the cell faces for the outward-pointing unit normal n.It is not difficult to show that there is a unique solution to the boundary-value problem (15) with vanishing mean, and we may define the apparent periodic conductivity It might also be shown that this conductivity satisfies the bounds and converges to the effective conductivity lim

Neumann boundary conditions
A third rather popular boundary condition is based on Neumann boundary conditions.Instead of the effective conductivity, the effective resistivity, the inverse of the effective conductivity, is approximated.Therefore, we fix a macroscopic heat-flux vector q, and seek the Neumann corrector  N,L ∶ Y L → R solving the boundary-value problem There is a unique solution with vanishing mean to this system of equations.Then, the apparent Neumann conductivity is defined via The key to the Neumann case is averaging the temperature gradient, as the boundary condition (19) fixes the mean value of the heat flux In the Neumann case, the defined apparent conductivity (20) satisfies the bounds and recovers the effective conductivity lim in the infinite-volume limit.
In general, the answer to the question which boundary condition to use depends on a number of factors.Most critical, to be emphasized again, is the fact that the defined apparent conductivities

2.3
A unified formulation of cell problems in homogenization

Generalities
There are some obvious similarities and differences between the three different apparent properties considered in the previous Section 2.2.The purpose of this section is to bring the problems into a common form to make them amenable to a unified solution strategy.For this purpose, it is convenient to recast the individual problems into variational formulations on suitable Hilbert spaces.To simplify notation, we will drop the explicit dependence of the correctors and the conductivities on the size L of the considered unit cell Y L .

Dirichlet boundary conditions
Let us start with the Dirichlet problem (11).We may express the Dirichlet corrector D as the (unique) minimizer of the variational problem where H 1 0 (Y L ) denotes the first-order Sobolev space of square-integrable functions which vanish on the boundary Y L of the cell.Introducing the weak gradient operator where the subscript D stands for "Dirichlet," and its formal adjoint mapping to the continuous dual space of H 1 0 (Y L ) and connected to the gradient D D by the equation we may express the Euler-Lagrange equation of the variational problem (25) in the form The difference between the Equation (25) and the original boundary value problem (11) is that the just derived equation ( 25) is a well-defined equation on a Banach (even Hilbert) space.The boundary conditions were absorbed into the definition of the spaces and of the operators D D as well as D * D .

Periodic boundary conditions
Let us turn our attention to the periodic problem (15).Its solution may be expressed via the quadratic optimization problem where H 1 per (Y L ) denotes the H 1 -closure of smooth functions which are periodic, together with their first derivatives, on the cell Y L and have vanishing mean value.
The similarities between the variational formulations ( 25) and ( 30) are rather striking.In fact, the formulations differ only by the considered ansatz space.Similarly, as we did for the Dirichlet problem (25), we may introduce the (weak) gradient operator and its formal adjoint which is defined via the prescription Although the definitions ( 26) and ( 31) look rather similar, the involved linear operators are different, essentially because they have a different domain of definition.In turn, their adjoints ( 27) and ( 32) are different, as well.
In analogy to the Dirichlet case, the first-order necessary condition of the variational problem (30) attains the form

Neumann boundary conditions
Last but not least, we turn our attention to the case of Neumann boundary conditions (19).The solution N arises as the minimizer of the variational problem The space H 1 # (Y L ) denotes the space of H 1 -Sobolev functions with vanishing mean.In contrast to the previous variational formulations (25) and (30), the integrand takes a slightly different form, owing to the prescribed heat flux.The problem (35) could be brought into a rather similar form utilizing the equivalent minimization problem by completing the square.Introducing the (weak) gradient operator for the Neumann case and its formal adjoint which is defined via the rule we may express the critical points of the Neumann problem in the form

Unified formulation
In the remainder of this section, we will introduce a mathematical abstraction layer which encompasses the continuous formulations of the cell problems for all three considered boundary conditions.Moreover, it will turn out to cover typical discretization schemes for these problems, as well, permitting us to carry out a unified analysis, and, subsequently, a unified design of numerical methods.Let Hilbert spaces H and V be given, together with a bounded linear operator s.t. the symmetric bilinear form defines an (equivalent) inner product on H, where (⋅, ⋅) V stands for the inner product on V.Moreover, suppose that a self-adjoint bounded linear operator is given which satisfies the bounds for some fixed positive constants  ± .Last but not least, suppose an element q ∈ V is given.Then, we consider the problem involving the adjoint D * ∈ L(V, H ′ ), defined via Existence and uniqueness of solutions to the abstract problem (45) are readily established.In fact, by definition of the adjoint operator (46), the problem ( 45) is equivalent to The left-hand side of Equation ( 47) defines a symmetric and coercive bilinear form on the Hilbert space H (due to the properties ( 42) and ( 44)), whereas the right-hand side gives rise to a bounded linear form.As a consequence, the existence and uniqueness of solutions is a consequence of Riesz' representation theorem for Hilbert spaces. 67Before turning our attention to numerical strategies for solving the problem (45), comprising Section 3, let us discuss how the different boundary conditions fit into the scheme.An overview of the different realizations of the boundary conditions within the described framework is given in Table 1.Establishing the defining properties for the boundary conditions at hand is straightforward except for the coercivity statement required for the bilinear form (42).The coercivity follows from suitable Poincaré inequalities.More precisely, (45).

Boundary condition
there is a positive constant C, depending only on the cell Y L , s.t.
as well as holds, see Alt. 67 The Poincaré inequality (48) implies coercivity for Dirichlet boundary conditions, whereas the second Poincaré-type inequality (49) may be used to establish coercivity for periodic and Neumann boundary conditions.
We will see in the subsequent Section 3 that the unified description introduced in the section at hand also encompasses discretized cell problems and different boundary conditions, at least provided a suitable discretization is used.In particular, the analysis extends to these cases as well, providing a unified treatment of the continuous as well as discrete problems irrespective of the boundary conditions.

A unified formulation of Green's operator methods
In this section we revisit the problem (45) that we wish to solve for  ∈ V. To set up a computational scheme, we utilize the assumption that the bilinear form ( 42) defines an inner product on the Hilbert space H. Riesz' representation theorem 67 asserts the existence of a bounded linear operator G ∈ L(H ′ , H), s.t. the identity The definition of the inner product (51) implies the identity Moreover, by definition of the adjoint operator ( 46), the Equation ( 52) is equivalent to As both f ∈ H ′ and  ∈ H are arbitrary, we thus deduce the identity In words, the operator G is the inverse of the (negative) "Laplacian" D * D, and we refer to this operator as Green's operator.
We thus observe that for any parameter  0 > 0, a field  ∈ H solves the problem (50) precisely if it is a solution to the Lippmann-Schwinger equation 68 which we may also express more compactly in the form with the notations G 0 =  −1 0 G and  0 =  0 Id.In engineering practice, the equivalent formulation of the Lippmann-Schwinger equation ( 57) on the temperature gradient  = D, is encountered rather frequently.However, the temperature-based formulation (57) permits implementations with a lower memory footprint.The Lippmann-Schwinger equation ( 57) is a fixed-point equation.For any reference conductivity  0 , the Lippmann-Schwinger equation ( 57) has a unique solution  * ∈ H, which is also a solution to the Equation ( 50) we started out with.The fixed-point form (57) naturally induces an iterative scheme with arbitrary initial value  0 .The iterative scheme ( 59) was introduced by Moulinec-Suquet 4 in the context of periodic cell problems in small-strain quasi-static mechanics and was subsequently analyzed by Eyre-Milton 69 for conduction problems.With similar arguments, the more general iterative scheme (59) may be analyzed to provide the error estimate with the prefactor which involves the upper and lower bounds (44) on the conductivity tensor.The error estimate ( 60) is formulated in the norm corresponding to the bilinear form (51), which is actually the natural norm to measure convergence for this kind of problem.We observe from the estimate (60) that the iterative scheme (59) converges provided the contraction factor ( 61) is strictly less than unity.This is the case provided the reference conductivity is chosen according to the rule  >  + ∕2.
The contraction factor ( 61) is minimized for The corresponding contraction estimate (60) becomes Thus, we observe that this contraction estimate involves solely the constitutive parameters, entering indirectly via the bounds  ± , and the initial guess  0 .Typically,  0 ≡ 0 is chosen, s.t.this influence factor does not enter, as well.We wish to emphasize that in this general setup, the (bounds on the) convergence rate is actually independent of the spaces H and V as well as the gradient operator D. Of course, these items enter the "measurement device," that is, the norm.Still, returning our focus to different boundary conditions, we observe that each of these boundary conditions gives rise to the same convergence rate when using these Green's operator methods.A discussion of a numerical method is not complete unless a proper convergence criterion is specified.Based on the identity we notice the representation of the natural residual associated to the governing Equation (50).Thus, we are led to the convergence criterion There is a number 30,31,33,70 of schemes similar to the basic scheme ( 59), but with an improved convergence behavior.Classically, these methods utilize the FFT, that is, the evaluation of the involved Green's operator can be considered exact (up to machine precision).However, in our study of computational homogenization methods based on Green's operators, we will also allow for less accurately evaluated Green's operators.Some accelerated numerical methods hinted at earlier are more robust to such an inaccurate evaluation.As a representative of this class, we will consider the Barzilai-Borwein (BB) method, 29 originally introduced as a quasi-Newton method.In our context, the BB method leads to the iterative scheme with a reference conductivity  k that changes from iteration to iteration according to the rule which amounts to a two-point approximation of the BFGS quasi-Newton scheme. 29he pseudocode for the implemented basic scheme with Barzilai-Borwein step size control is given in Algorithm 1.

The explicit jump discretization
Various discretization schemes were proposed to solve the Lippmann-Schwinger equation on a discrete grid.For reasons of compatibility with the discrete finite difference operators in the TT format-which we will discuss in more detail in the following section-an implementation using finite difference operators is advantageous.When solving the Lippmann-Schwinger equations associated to a general finite difference scheme, solution artifacts like checkerboards might occur, leading to nonphysical local fields.Additionally, general finite difference schemes may introduce an asymmetry of the solution field, which may impair enforcing the boundary conditions properly.To avoid these difficulties, we employ the EJ discretization, 27,71 which, by a change of variables, can be treated like a finite difference discretization with a modified evaluation of the material law.To simplify the presentation, we restrict to a three-dimensional cubic domain Y = [0, L] 3 which we discretize into N 3 equally-sized voxels of edge length h = L∕N.In fact, the extension to general rectangular domains is possible at the expense of more involved notation.The EJ discretization employs this voxel grid with temperature values living in the center of the voxel and collected in a vector  ∈ R N  .Heat fluxes q ∈ R N  live on the interfaces between voxels.The evaluation of the finite differences on boundary voxels differs between the boundary types, but interior voxels are treated in the same way for all boundary condition types.Figure 1 shows a sketch of a single voxel for the EJ discretization.For exposition, the grid is shown in two dimensions.The purpose of the EJ discretization is solving the balance of the heat flux for every voxel on the discrete grid.This is equivalent to the condition that the sum of all heat fluxes on each voxel's interface vanish.Like the heat flux values, the discrete (periodic) temperature gradients  ∈ R N  ,  = D h  are located on the faces of the voxel adjacent to the voxel needed to compute the finite difference gradient D h .The components of the discrete temperature gradient  on the interior are computed via For example, for the x-direction the gradient  x o,p,q is located halfway on the line connecting the locations of  o,p,q and  i−1,j,k on their respective voxel centers.Discrete heat fluxes  =  with a conductivity tensor  = diag(, , ) are located on the voxel faces at the same position as the respective temperature gradients  and are computed via a harmonically averaged material law defined on each of the voxel centers.To finish defining the discretization scheme (67), we introduce the divergence of the heat flux D * h  .For the EJ discretization, we compute the divergence of the heat flux as where the values obtained from an application of the divergence operator D * h  live on the same locations as the temperature values .
Defining the differential operators for the gradient (70) and the divergence (72) in the presented manner makes their evaluation algorithmically equivalent to a finite difference discretization, and the only difference lies in the harmonically averaged material law (71).Most crucially, the computation of the Laplacian Δ h = D * h D h is identical to that of the conventional eight-point finite difference Laplacian in 3D, which we are going to use for computing Green's operator G h = Δ −1 h in the next section.
In the discrete unified framework with Green's operator the application of the boundary conditions on the domain edges is handled by the differential operators.This gives rise to different discrete differential operators for each type of boundary condition.The algorithmically simplest case comprises periodic boundary condition with the gradient operator D P,h , where finite differences on the domain edges are handled like voxels on the interior, but the periodically neighboring voxel is chosen.Therefore, for each direction, the gradient field has the same number of entries as the temperature field N  , with the full gradient field having the threefold number of entries N  .An overview of the number of grid points for all three BC types in given in Table 2.With these tools at hand, we formulate the solution scheme for the periodic case as the discrete analogue of the variational problem ( 30) TA B L E 2 Domain sizes for discrete temperature and gradient fields for different boundary conditions and the EJ discretization.

No. of temperature points No. of gradient points
Periodic with Green's operator and where  refers to the prescribed mean temperature gradient constant on all voxels.For Dirichlet boundary conditions, Figure 2A shows a two-dimensional sketch of the grid where the boundary values for the temperature are prescribed by means of so-called "ghost points" which are denoted in brackets.Therefore, we compute the gradients D D,h on all faces of voxels belonging to the interior of the domain.Consequently, for a fixed number of temperature points, we need to compute one more gradient value in direction of the gradient evaluation, see Table 2.By setting the value of the ghost points, we set the value of the Dirichlet boundary on the domain edge, although we restrict to zero boundary conditions.
The ghost points enable us to compute D * hD h (⋅) on any discretized point which amounts to the computation of the 2nd order Laplacian with Dirichlet boundary conditions, a property we will use in the next section to compute Green's operator.In contrast to other approaches, 14 we do not need to introduce internal forces to apply the boundary condition.The discrete solution scheme with an applied mean temperature gradient  for the Dirichlet case is with Green's operator For Neumann boundary conditions, Figure 2B shows a two-dimensional sketch of the domain for a small problem.By setting the values of the gradients on the boundaries, we implement the Neumann boundary condition.For a fixed number of temperature points, Neumann boundary conditions actually need one gradient value less in direction of the gradient evaluation than periodic boundary conditions, as the gradients on boundary are fixed.The discrete solution scheme differs from the schemes for periodic or Dirichlet BCs.As prescribing zero heat flux on the boundary in combination with a constant mean temperature gradient allows for the trivial solution  = 0, in the case of Neumann boundary conditions, we prescribe a mean heat flux.Therefore, the discrete solution scheme for Neumann boundary conditions becomes with Green's operator and the mean heat flux q.Thanks to the unified framework, the solution schemes for all three BC types, that is, periodic (75), Dirichlet (77), and Neumann (79) BCs, are similar even in the discrete case, save for the differential operators and the type of the applied mean field.The established EJ discretization (73) fits into the abstract scheme of Section 2.3, as the operators D and D * are adjoint by construction and the non-degeneracy condition 67(section 12) holds as a consequence of a discrete Poincaré inequality.

3.3
Computing Green's operator with Fourier methods

Generalities
Applying Green's operator is a bottleneck for the runtime when solving the cell problem via the Lippmann-Schwinger approach (57).Taking a step back to the continuous case, we recall the (continuous) Lippmann-Schwinger equation which we rewrite as

Periodic boundary conditions
We expand any sufficiently smooth and Y -periodic function  into a Fourier series in terms of the imaginary unit i, that is, i 2 = −1, with Fourier coefficients Then, applying the periodic continuous Green's operator G P = ( D * P D P ) −1 from Equation ( 82) amounts to a multiplication of the Fourier coefficients by the inverse eigenvalues of the Laplace operator. 60n the discrete case, we employ an analogous approach. 60We use the FFT to compute the (discrete) Fourier transform of the right-hand-side of the discrete periodic Lippmann-Schwinger equation on the EJ grid, 27 that is, When using such finite differences instead of derivatives, the action of the Green's operator from Equation (86) in the periodic case 17,60 becomes Thus, periodic boundary conditions are covered naturally via the Fourier series and, thus, the discrete Fourier transform.

Dirichlet boundary conditions
To transfer this straightforward approach to Dirichlet and Neumann boundary conditions, we employ appropriate sine and cosine transforms instead.These transforms are real-valued relatives of the complex-valued Fourier transform, with either only the sine or cosine terms as their basis functions.The sine transform intrinsically encodes Dirichlet boundary conditions, whereas the cosine transform enforces zero Neumann boundaries.The continuous sine series for a smooth and real function  takes the form with the sine coefficients We express Green's operator with Dirichlet boundary conditions Subsequently, we will use the discrete sine transform DST-I in all three dimensions as the discrete analog to the sine series. 59,60We solve the discrete Lippmann-Schwinger equation for Dirichlet boundary conditions on the finite difference grid via the coefficients of the DST-I )) 2 . (92)

Neumann boundary conditions
Last but not least, we take a look at Neumann boundary conditions which we enforce using cosine series/the cosine transform.Expanding a function  into a cosine series with cosine coefficients permits enforcing Neumann boundary conditions.
In this case, the application of Green's operator reads Then, in analogy to the periodic and Dirichlet case, we use the discrete cosine transform DCT-I to compute the action of Green's operator for Neumann boundary conditions on the EJ grid via scaling the cosine coefficients To conclude, Equations ( 92), (87), and (97) permit to compute the update step for the FFT-based method of applying Green's operator for Dirichlet, periodic, and Neumann boundary conditions, respectively.In the next subsection, we will show an alternative strategy using the low-rank TT-based approach which requires no Fourier-like transforms.

Computing Green's operator on a tensor train
As the operators D and D * for finite differences (and therefore for the EJ discretization) are linear mappings on finite-dimensional spaces, they may be recast in the form of matrices.Then, by definition ( 52) Green's operator may also be represented as a matrix.To implement the iteration scheme (67), we may thus apply Green's operator to the right-hand side of iteration scheme via a matrix-vector product.As Green's operator is fixed up to a scaling constant in our solution scheme, see Algorithm 1, Green's operator may be precomputed and cached.Green's operator on an N 3 grid stored as a full matrix requires N 6 entries in memory-therefore, the outlined approach is virtually unfeasible using a conventional full matrix format.Moreover, the matrix is fully populated, and sparse matrix techniques do not apply.Therefore, we use the QTT format 55 as a low-rank approximation method to reduce the memory footprint of the matrix.The QTT format-the application of the TT format to low-dimensional problems-allows for an efficient storage of certain matrices like the multidimensional finite-difference Laplacian and its inverse.For a p-dimensional array-like structure T ∈ R n 1 ×n 2 ×•••×n p with n k entries in the kth dimension, the TT representation approximates a given p-dimensional array as a product of p three-dimensional arrays, which are called TT cores G i .The size and memory footprint of these tensors is governed by the individual cores' TT ranks r i .For conventional arrays, exactly n p array entries need to be stored, provided n 1 = n 2 = • • • = n p = n holds.Low TT ranks lead to a significant memory reduction by storing only the TT cores.In case of the QTT format and for fixed TT rank, the memory required grows only logarithmically in the number of entries of the original array.The QTT format enables a low-rank representation of linear operators, provided their structure allows for such a representation, as for low TT ranks, the memory requirements are very low.Kazeev and Khoromskij 43 derived an explicit expression for the inverse finite difference Laplacian in one spatial dimension with minimum QTT ranks.For the higher-dimensional Laplacian, no explicit representation of the inverse is not known, but an approximation of the inverse high-dimensional Laplacians with low TT ranks is possible.We construct higher-dimensional Laplacians as an outer product of one-dimensional Laplacians using a so-called Kronecker product in the QTT format, 43 making it a viable to compute Laplacians for non-cubic rectangular domains as long as the numbers of voxels per edge length are powers of two.For simplicity, in this contribution we restrict to cubic domains.The boundary conditions, which require using completely different Fourier transforms in the approach previously discussed, are encoded directly in the one-dimensional Laplacians for this QTT-based approach. 43Nonetheless, some care has to be taken concerning the boundary conditions.Apparently, a matrix needs to be regular to permit an inverse.For periodic and Neumann boundary conditions, the Laplacian matrix has a non-trivial kernel consisting of constant fields, that is, multiples of the field 1.
To handle these cases, we make use of the following construction.Suppose a symmetric and positive semidefinite matrix A ∈ R n×n is given with a one-dimensional kernel represented by scalar multiples of the non-zero vector v.Then, for any positive scalar c > 0, the shifted matrix is positive definite, and thus invertible.Suppose we wish to solve the equation for a right-hand side b that is orthogonal to the kernel of A, that is, satisfies v T b = 0.Then, we claim that the Equation (100) is solved by x = Ã −1 b.In fact, by definition, the vector x satisfies Multiplying this equation by the vector v, we observe where we used the symmetry of A and the fact that v lies in the kernel of Thus, in view of the representation (101), the vector x satisfies the Equation (100), as desired.
For the problem at hand, periodic and Neumann boundary conditions, the matrix A corresponds to the (negative) finite-difference Laplacian and the kernel is given by constant fields, that is, we v ≡ 1.As the shifting strategy (99) preserves an existing low-rank structure of the matrix A, we may thus obtain Green's operator in these cases from the shifted Laplacian (99) where the constant c is chosen to lie within the non-zero eigenvalue spectrum of the operator A.
Using methods like the LU decomposition 72 to compute the inverse of the Laplacian stored in the QTT format is not suitable for computing Green's operator, because accessing individual entries of QTT matrices is rather expensive.Thus, to compute an inverse of a matrix in the QTT format efficiently, it is imperative to rely on a strategy that does not require access to matrix entries.One such method is the computation of inverses via exponential sums.In the spirit of the work by Beylkin et al., 73 this approach, which was applied to low-rank approximations 74,75 and especially the QTT format 56 before, inverts a positive definite matrix A by a quadrature of the integral representation of the matrix inverse, involving a positive parameter a.In the discrete case, we use a method by Khoromskij and Oseledets, 56 called the t k -quadrature, which approximates the inverse via and involves a parameter M. To evaluate this sum, the matrix exponential exp(−t k A) needs to be computed.Following Khoromskij and Oseledets, 56 we employ the scaling and squaring 56,76 approach.We select a sensible m, so that e A∕m in is simpler to compute than the original exponential.Then, we compute e A∕m by a Padé approximation.Matrix multiplications and matrix-vector products, which are unavoidable during the inversion of high-dimensional Laplacians using the scaling and squaring approach, increase the ranks of QTT matrices of the resulting matrix compared to the factors.We perform a so-called TT recompression, 15 where we prescribe a maximum rank to prevent excessive memory use.The computation of the matrix exponential of a given matrix in the QTT format is summarized in Algorithm 2.
Algorithm 2. QTT matrix exponential using the scaling and squaring method 56,76 Require: QTT matrix A, QTT identity I, Taylor cutoff M exp Ensure: QTT exponential matrix e A s ← max ( 0, ⌈log 2 ‖A‖⌉ ) k CA + I end for for l = 1, … , s do ⊳ Squaring C ← CC TT compression with predefined rank limit r max,G end for return C Using this matrix exponential, we eventually compute the inverse of the Laplacian using the truncated series (104).

Setup
In this section, we investigate and compare the numerical performance of applying Green's operator using both the TT and the Fourier approach.To this end, we first investigate the numerical behavior of the TT inversion algorithm and then turn our focus to the differences in performance for the different boundary conditions.Using the Green's operators for different boundary conditions, we will investigate the effect of different boundary conditions on the local solution fields of the cell problems.Finally, we will investigate how different boundary conditions interplay with different types of microstructure ensembles.To apply Green's operator on the EJ grid, the discrete field  needs to be transformed by a (real) discrete Fourier transform for which we employ the FFTW 77 library.We use a DST-I in all three dimensions as the discrete sine transform for Dirichlet boundary conditions.In analogy, we use the DCT-I as the discrete cosine transform for Neumann boundary conditions.For periodic boundary conditions, we employ the three-dimensional FFT.The experiments were performed using our own in-house tensor train library and homogenization solver for heat conduction, both written in Python with Cython extensions for parallel execution of performance-critical sections of code.Runtimes were recorded on an 8-core Intel i7-10700 CPU with 64GB of RAM.

Numerics of Green's operator
Our solution scheme (67) for the homogenization problem involves Green's operator serving as a preconditioner for the solution scheme.A sensible accuracy of the preconditioner is therefore crucial for fast convergence.The accuracy of the preconditioner applied using the Fourier approach from Section 3.3 is only limited by floating point accuracy.The QTT-based approach, on the other hand, involves an approximation of the inverse Laplacian using the quadrature as well as the scaling and squaring scheme.To investigate the accuracy of the QTT-based Green's operator with different boundary conditions, we compute the inverse of the Laplacian with the respective boundary conditions in the QTT format.To start, the Laplacians have rather low maximum QTT ranks of five for periodic and Dirichlet boundary conditions, and seven for Neumann boundary conditions.Then, to obtain the inverse of the Laplacian, that is, Green's operator, we use the quadrature (104) with the matrix exponential computed via the scaling and squaring approach, see Algorithm 2. We use an assembly-free matrix-vector product 78 for applying Green's operator stored in the QTT format and store the right-hand side of Equation ( 67) as a full vector.To evaluate the accuracy of computing the approximated Green's operator G TT , we use the error Error err TT of Green's operator obtained by inversion of the Laplacian using the QTT approach on a 128 3 EJ grid for three types of boundary condition and for multiple limits on the TT ranks.Here, r max is the rank limit enforced after inversion.
giving us a measure of the accuracy of the inversion, involving the identity matrix I TT and the Laplacian Δ TT in the QTT format.Various parameters influence the accuracy of the inverse.For the cutoff of the Taylor series (see Algorithm 2), we used a value of M exp = 15. 76We selected a high value for the parameter M quad = 30 as this gives the most accurate results for all cases.Reducing this value does not influence runtimes significantly because the inversion is only performed once per computation run. 43s matrix products of two objects in the TT format multiplies their TT ranks, a sustained TT compression is needed after every product to keep the TT ranks from growing in an uncontrollable way.Therefore, the limit on the TT ranks r max of the matrices during evaluation of the matrix exponential and during evaluation of the inverse represents another parameter that has an influence on the accuracy.Setting these parameters too low limits the accuracy but choosing the rank limit before a matrix product too high hampers performance because of the multiplicative nature of the ranks.To quantify the accuracy of the inversion, we evaluate Green's operator for the EJ discretization on a 128 3 grid for all three types of boundary conditions.In Figure 3, the accuracy of Green's operator for various TT rank limit imposed during and after inversion is shown.
For low QTT rank limits, we observe a comparatively high error for all boundary conditions, where the Neumann boundary condition shows the highest error.Increasing the rank limit reduces the error for all boundary conditions up to the rank limit r max = 32, where the periodic Green's operator reaches its lowest error.Green's operator with Neumann boundary conditions attains a similar level, but only for the doubled rank limit of r max = 64.Dirichlet boundary conditions on Green's operator show a larger error at a rank limit r max = 32, but at higher rank limits, solely the Dirichlet boundary conditions allow for a further reduction of the error.Therefore, on a 128 3 grid, a rank limit r max = 64 allows for a similar accuracy for all boundary types.Only the Laplacian with Dirichlet boundary conditions yields a positive definite matrix.Even using the slight perturbation does not allow for an inversion of the periodic or Neumann Laplacian with an accuracy equivalent to the Dirichlet case.
The evaluation time for a product of a QTT matrix and a full vector is affected by the ranks of the QTT matrix.Thus, the higher ranks needed to obtain an adequately low error lead to higher evaluation times for Neumann and Dirichlet boundary conditions.In Figure 4, runtimes of the matrix-vector product of the QTT Green operator are plotted over the rank limit.We observe a strong influence of the rank limit r max on the evaluation times for all three types of boundary conditions.The runtime to evaluate Green's operator remains steady for all types of boundary conditions, as we only enforce hard rank limits after the evaluation of Green's operator.Unfortunately, this also means that for larger domains, where we have to scale QTT rank limits with the edge length N for a consistent approximation error, the runtimes will increase too.Noticeably, for larger domain sizes, the evaluation of Green's operator with periodic boundary conditions is almost an order of magnitude faster than with Neumann boundary conditions.This effect is caused by the lower TT ranks needed to store the periodic Green's operator.For all types of boundary conditions, the runtimes level off for higher Runtimes of a single product of a QTT Green operator with a full vector on a EJ grid of size 128 3 .The rank limit of Green's operator is denoted by r max .

F I G U R E 5
Evaluation times of a single application of Green's operator on an EJ grid of size N 3 for the QTT-and FFT-based approaches and all three types of boundary conditions, respectively.rank limits.This is caused by the TT recompression, keeping the ranks as low as needed, as long as they remain below the rank limit.For instance, in the periodic case, the runtimes barely increase for rank limits r max larger than 64.In this case, the actual ranks of the TT cores are below the rank limit and do not have detrimental effect on the evaluation times.For Neumann boundary conditions, the actual ranks of the Green's operator are higher.Consequently, an increase of the rank limit above 64 still leads to an increase in runtime.For reference, in Figure 5, we show the runtimes for a single application of Green's operator using the Fourier approach for all three types of boundary condition, where the evaluation time includes the Fourier transforms of the field into Fourier space and back into real space.We observed that throughout, the runtimes of the Fourier-based approach are lower by two orders of magnitude compared to the QTT-based approach.The evaluation times surpass more than 100s for even medium-sized domains like N = 256.For this domain size, the evaluation for the FFT-based Green's operator are still below 1 s.A reason for this speed-up is that the efficient parallelization of the Fourier transform has been the subject of extensive work 77,79 and the application of Green's operator (82) is parallelized without much effort.However, parallelizing operations on TT matrices is not as straightforward due to the memory structure and interconnected nature of the TT cores, especially the very small TT cores in QTT, see Daas et al. 80 In addition, in Figure 5 we observe the that the runtimes of the FFT-based approach scale linearly with the domain edge length, similar to the QTT-based approach.

F I G U R E 6
Comparison of the number of iterations to a residual below tolerance 10 −5 for different domain resolutions N and types of boundary conditions between the QTT-and FFT-based Green's operators.For the QTT-based Green's operator, a rank limit r max = 128 was selected.
Next, we investigate the convergence of the solution scheme using the Barzilai-Borwein step size control (Algorithm 1) where Green's operator serves as the preconditioner.We consider the two approaches of applying Green's operator and compare all three types of boundary conditions.To this end, in Figure 6, the number of iterations needed to reach a residual below 10 −5 over the domain resolution N 3 is shown.The microstructure is a single spherical inclusion, sketched in Figure 9, filling a volume fraction of 0.15 the cell and centered in the cell with a medium material ratio  1 ∕ 2 = 10.
We observe that, in general, the domain size does influence the iteration count for the QTT-based Green's operator.For the QTT-based approach this means that increasing the domain size reduces the quality of the preconditioner.The number of iterations for the QTT-based method exceeds those of the FFT-based method for the same boundary condition by a factor of two, at least, and increases linearly with the domain size N.The reduced accuracy of the preconditioner is caused by the fixed TT rank limit r max = 128.If we do not fix the rank limit, evaluation times for domains with resolutions N larger than 64 become unsustainably large.As discussed above, the QTT-based inversion scheme used to obtain Green's operator has only limited accuracy.In our case, the accuracy of the preconditioner is still sufficient to ensure convergence for all our test cases, but fails to generate minimum iteration counts. 72The TT Green's operator with Neumann boundary conditions shows distinctly higher iteration counts than periodic and Dirichlet boundary condition in spite of this operator reaching the same accuracy as the periodic operator during the inversion procedure (Figure 3).Instead, the iteration count for the periodic and Dirichlet case differ only by a single iteration, even though the inversion error of the Dirichlet Green's operator is almost two orders of magnitude lower.
For the FFT-based preconditioners, the boundary conditions have only a minor influence on the respective iteration counts.In addition, the iteration counts for the FFT-based preconditioners are considerably lower.Both effects, the influence of the boundary conditions and lower iteration counts, are caused by the exact inversion of the Laplacian when using the FFT-based approach.At last, the number of iterations does not increase when the domain size N is increased which also is an effect of the decreased accuracy of the inversion of the preconditioner for the QTT-based approach.Therefore, Green's operator for the FFT-based method is an exact preconditioner for the Lippmann-Schwinger equation.In the Moulinec-Suquet method, it is sometimes observed that the effective properties converge but the solution field does not.To rule out such a phenomenon, we investigate the residuals for our solution method.The residual at hand (Algorithm 1) quantifies convergence of the local solution field.For exposure, we collected the residuals in Figure 7 for a domain size of N = 64.We observe that for all approaches, the residual declines in the manner typical for the Barzilai-Borwein method 30 and eventually falls below the selected tolerance.

F I G U R E 7
Comparison of the residuals for different types of boundary conditions for a resolution N = 64 between the QTT-and FFT-based Green's operators.

F I G U R E 8
Comparison of the evaluation time up to a residual below tolerance 10 −5 for different domain resolutions N and types of boundary conditions between the QTT-and FFT-based Green's operators.For the QTT-based Green's operator, a rank limit r max = 128 was selected.
Finally, we compare the runtimes for the evaluation of the whole solution scheme in Figure 8.The detrimental effect of the slower application of the QTT-based preconditioner is reinforced by the higher iteration count needed which in turn is caused by the lower accuracy of the preconditioner.For both the QTT-based and the FFT-based preconditioner, the evaluation time in the Neumann case is highest.This is caused by the slightly higher number of grid points that need to evaluated, as the higher evaluation time is also present for the FFT-based preconditioner.Throughout, the FFT-based method is two orders of magnitude faster.For instance, finding the converged solution on a 128 3 grid with periodic boundary conditions took 14.5 s using the FFT-based preconditioner, but 1252 s using the QTT-based preconditioner.
In the following subsections, we will investigate the effect of boundary conditions on homogenized quantities.Given the relative slow evaluation of the QTT-based preconditioner, we will employ the FFT-based approach and discard the QTT-based approach.
To conclude, we observe that to precondition a symmetric and positive definite problem, the preconditioner should have a similar structure as the system matrix itself.We can ensure this by employing the same discretization but with constant coefficients for the preconditioner (which is the EJ discretization in this work).Here we observe a similar effect when only approximating the preconditioner, that is, evaluating the preconditioner accurately turns out to be critical for optimum numerical performance.
F I G U R E 9 Sketch of a section a single spherical inclusion (left) and a section through Hashin's coated sphere (right).

Homogenization with boundary conditions
Only for sufficiently large microstructure cells, the influence of the boundary conditions on the apparent properties becomes negligible. 2 Using our unified Green's operator based approach, we evaluate how the boundary conditions influence the apparent properties of the material and study the approximation capabilities of the effective properties.
To investigate the quality of approximation of the homogenization scheme with different boundary without resorting to the limit case, we use Hashin's coated sphere, a material with a known analytical solution for the cell problem. 81In the case of thermal conductivity-as sketched in Figure 9-Hashin's coated sphere consists of an interior sphere with conductivity  1 and radius r 1 that is expanded by a coating of conductivity  2 , so that the thickened sphere has a radius r 2 .The coated sphere is embedded in a material with a conductivity  eff .The values of the conductivities are computed from a prescribed effective conductivity.Given a non-negative material ratio , we set for the conductivity of the sphere.Using the abbreviations f 1 = (r 1 ∕r 2 ) 3 and f 2 = 1 − f 1 , we select the conductivity of the coating as where the effective conductivity of the whole domain is equal to the conductivity of the embedding material  eff .In other words, by choosing the parameters deliberately, the coated sphere becomes "invisible" in view of the effective properties.
To validate the results of our implementation with analytical results, in our numerical homogenization scheme, we prescribe the effective conductivity of the embedding  eff an and compute the effective conductivity  eff num .The relative error  rel then computes as As each coordinate direction shows the same behavior because of the isotropy of the effective conductivity, we will only evaluate the values in a single direction.The microstructure with a single Hashin's sphere is perfectly periodic and has a deterministic distribution of the microstructure.Therefore, we carry out a resolution study by performing the periodic homogenization on the domain Y = [0, L] 3 of fixed size L = 1 m, with each dimension discretized into a varying number N of voxels.In Figure 10, the relative error is plotted for all three edge boundary types for various resolution steps N.
A low material ratio  of 10 −2 , corresponding to a sphere of low conductivity coated by material of high conductivity, provides a relative error for the periodic boundary conditions that is lower by almost an order of magnitude when comparing to the other two boundary condition types.In this case the relative error scales proportionally to the voxel size h = L∕N.This effect has been called superconvergence, as the relative error of the local fields only converges with the square root of the voxel size h 1∕2 . 82Nonperiodic boundary conditions can disturb the homogeneity of the solution field. 9

F I G U R E 10
Relative error (107) for two different values of the material ratio  and for radii 2r 1 = r 2 = 0.4 for Hashin's sphere.
In case of Hashin's sphere, the field surrounding the sphere is perfectly homogeneous.Therefore, a better rate of convergence of the periodic boundary conditions is plausible, even though in the continuous case, all boundary conditions lead to the same solution.
Interestingly, the convergence behavior of the relative error for Neumann and Dirichlet boundary conditions is virtually indistinguishable, yet an order of magnitude larger than for periodic boundary conditions.For a material ratio  of 10 2 and periodic boundary conditions, the relative error is reduced by an order of magnitude for low resolutions N compared to both Neumann and Dirichlet boundary conditions.In contrast, the change in material ratio barely affects the relative errors for Neumann and Dirichlet boundary conditions.Whereas for periodic boundary conditions the proportionality of the error to the voxel size h is lost for smaller grid sizes, Neumann and Dirichlet still scale proportional to h, though on a higher level.
In Figure 11, the iteration count for the solution schemes is shown for all three types of boundary conditions and at different resolutions.When compared to the single (uncoated) sphere in Figure 6, we notice that the considerably higher material ratios lead to higher iteration counts.Even though the iteration counts for the prior example with lower material contrast remained constant for an increase in resolution, the iteration count slightly increased for higher resolutions.This is the case for all types of boundary conditions.We observe that Neumann boundary conditions lead to the highest iteration counts but remain relatively constant, whereas periodic and Dirichlet boundary conditions lead to a more restricted rise in iteration count for higher resolutions.In contrast to the iteration counts for a single sphere in Figure 6, the iteration counts increase, even though iteration counts are bounded on a low level, independently of resolution with low resolution yielding even lower iteration counts.

Influence of microstructure and boundary conditions for stochastic materials
Obtaining periodic microstructure images is challenging unless the microstructures are generated.Three-dimensional images of real microstructures, for example, using micro-computed tomography (CT), generally do not have a periodic microstructure.Computing the apparent properties with periodic boundary conditions therefore introduces some error when features of the microstructure are cut off at the volume element's edges.Therefore, nonperiodic, that is, Dirichlet and Neumann, boundary conditions have the potential to lead to apparent properties that are closer to the effective properties than periodic boundary conditions.We consider spherical fillers with a conductivity 1.5 W∕(K m) where the spaces between the spheres inside a matrix material with a conductivity of 0.195 W∕(K m). Figure 12 shows different realizations of this microstructure type on volume elements with variable size.The smallest microstructure (L = 16 m) consists of a single whole sphere as well as three partially cut spheres.The largest microstructure (L = 512 m), on the other hand, contains around 122,500 spheres.We generate these microstructures 83,84 with a predetermined volume fraction of filler material.Then, we solve the homogenization problem in a unified manner with periodic, Neumann and Dirichlet boundary conditions for multiple types of ensemble configurations.
For each microstructure boundary type, we generate microstructures of varying edge length L with mono-sized spheres, so that the number of spheres in the domain increases proportional to L 3 .Then, we generate 20 realizations for each of these combinations to account for stochastic effects of the microstructure generation process.
We use the sequential addition and migration (SAM) 84 algorithm for generating the sphere packings.The SAM algorithm is able to generate packings of spherocylinders, in general, and spheres may be considered as a special type of spherocylinder.To generate periodized microstructures of non-overlapping spheres, potential overlaps with spheres in periodic repetitions of the domain need to be considered, as shown in the sketch in Figure 13A where the spheres belonging to the original domain are filled in a darker shade.This microstructure has spheres intersecting the edges of the domain.
Figure 14 shows the computed effective conductivities for a loading in x-direction for three smaller domain sizes L with a constant number of voxels per sphere and voxel size h = 1 m.The effective conductivity  eff of 0.2895W∕(K m) was computed from the periodic microstructure with periodic boundary conditions on the largest domain with a length L of 512 m.For the periodic microstructure (Figure 13A), the mean of the effective conductivity for a very small microstructure L = 16 m is very close to the effective properties of the reference solution with L = 512 m, underlining the fact that, for periodic microstructures, periodic boundary condition reproduce the effective properties best.In contrast, we observe that for Dirichlet boundary conditions, the effective conductivity is significantly overestimated.Moreover, the variation of the results is higher.This overestimation is rooted in the varying number of cut spheres, which have a larger influence for the Dirichlet case.The results for Neumann boundary conditions on the other hand are very close to the results of periodic boundary conditions even for small microstructures.In Figure 14, for L = 32 m the values of periodic and Neumann boundary conditions are virtually indistinguishable.
Subsequently, we investigate microstructures extracted from a larger, periodically generated domain.This strategy is intended to mimic results of image acquisition techniques from real-world materials.In Figure 13B, this procedure is depicted for a single direction, that is, we generate a periodized microstructure and extract only half of the generated data in x-direction.
For the study at hand, we perform the extraction procedure in all three space dimensions.The extraction leads to ensembles where the faces of the microstructures are nonperiodic, and spheres can cut these faces.For lower microstructure sizes up to 64 m, we observe in Figure 14 that the variance of the conductivity is larger for all types of extracted microstructures compared to the periodic microstructures.The increased variance is a direct consequence of the increased variation in the realized volume fraction caused by the extraction process.In fact, the periodically generated microstructures have a fixed volume fraction which is accurate to machine precision.In contrast, the extraction process does not account for the volume fraction, introducing randomness at this stage.Consequently, the effective conductivity is determined by the volume fraction to first order, 2 which is responsible for the observed high degree of randomness of the apparent conductivities.Even though this behavior disturbs the apparent conductivities from quickly reaching the effective conductivity, the considered procedure serves as a model for the real-life procedure of cutting the microstructure from a larger sample of the material.We notice that, for this microstructure type and for a small domain, Neumann boundary conditions show results closer to the effective properties calculated for L = 512 m than periodic boundary conditions, because the cut spheres only influence one side of the domain.As the spheres cut the faces, but do not have corresponding periodic images on the other side of the domain, the periodic boundary conditions introduce an error by creating sharp edges, leading to singularities in the solution fields.A different point of view considers the cut spheres as another species of particles which are introduced into the composite.Neumann boundary conditions do not couple two opposite faces and, for this type of ensemble, estimate the effective conductivity more accurately than periodic boundary conditions.For extracted microstructures, Dirichlet boundary conditions drastically overestimate the conductivity, because Dirichlet boundary conditions fix the temperature on the domain edges, therefore making the material appear more conductive than the considered microstructure actually is.
The last microstructure type with spherical inclusions under investigation are hard boundaries.In this case, spheres are not allowed to intersect the faces of the cells which is ensured in the microstructure-generation process via additional collision checks with the faces.The result is illustrated for a single direction in Figure 13C.Inspecting Figure 14, we observe that, on small domains, hard boundaries in all directions provide the lowest variance of the results compared to other microstructure types.Especially Dirichlet boundary conditions perform better than for the other considered cases, albeit they still overestimate the conductivity for all but the largest unit cells.Apparently, Dirichlet boundary conditions benefit from inhibiting spheres from intersecting the cell faces.For Neumann boundary conditions, the overestimation is stronger than for periodic microstructures, whereas the overestimation using Dirichlet is less grave than for other microstructure types such as the extracted or periodic microstructures.
For large domain sizes, for example, L = 512 m in Figure 14, the apparent conductivity values align for all microstructure types for Neumann and periodic boundary conditions.With larger domains, the influence of the domain boundary on the apparent values decreases.Only for Dirichlet boundary conditions, even for the large microstructures, the apparent conductivity lies about five percent higher.On the other hand, Neumann and periodic boundary conditions coincide up to four decimal places.All in all, Neumann boundary conditions overestimate the effective conductivity for all but some particular cases, confirming prior findings 2 that showed Dirichlet boundary conditions provide an overestimation compared to periodic and Neumann boundary conditions.

Influence of boundary conditions for industrial-scale realistic microstructures
After considering matrix-inclusion composites, we turn our attention to another class of microstructured materials of interest, namely porous materials.These feature a non-trivial volume fraction of material with a gaseous phase.With FFT-based computational homogenization schemes, such materials turn out to be rather computationally challenging, as the discrete systems involve a rather large material contrast.In turn, corresponding stochastic microstructures often decorrelate on much larger scales than materials with moderate contrast.In experiments, this phenomenon typically implies that the measurements cannot be conducted on a representative volume element.In fact, due to the huge correlation length, the measurements on samples of manageable size then critically depend on the boundary conditions applied during the test.As an example, we refer to the results of Ettemeyer et al. 85 who compared measured effective thermal conductivities of sand grain materials with computational homogenization results of periodic sand-grain microstructures (and periodic boundary conditions), and found noticeable deviations.In the section at hand, we will reconsider their problem and investigate the influence of the boundary conditions.
In Figure 15, a three-dimensional image of the generated periodic sand core microstructure is shown.The microstructure consists of two solid phases, the sand grains in gray and the binder in red as well as the voids in between the sand grains, which are assumed to have the thermal conductivity of air.
The numerical values of the conductivities are listed in Table 3 for two different types of sand grains. 85or six samples of a generated sand grain microstructure, we evaluate a unidirectional prescribed temperature gradient for each microstructure and for each of these, we evaluate three boundary types, Neumann, periodic, and Dirichlet.
Figure 16 shows a temperature gradient field for the three different boundary condition types.There are some subtle but notable differences in the local fields.In comparison to the periodic microstructure shown in Figure 16B, the Dirichlet boundary condition shown in Figure 16C leads to an increase of the temperature gradient in vicinity to the edges of the domain.This change underlines that the Dirichlet boundary condition acts like a highly conductive wall.Periodic boundary conditions, on the other hand, provide an undisturbed, that is, homogeneous, field on the domain boundaries.The temperature gradient field evaluated with Neumann boundary conditions shown in Figure 16A, shows high gradient values on the domain faces where the air phase of the sand core meets the boundaries that have a normal in the loaded x-direction.This effect is introduced by the fixed zero heat flux on the boundary that needs to compensate for the applied mean heat flux for Neumann boundary conditions.The temperature gradient field on the interior of the domain remains largely unaffected by the boundary conditions having a limited range of influence.In fact, this leads to a decreased influence of the boundary conditions with increasing size of the unit cell.Only for periodic boundary conditions, we observe a slightly higher temperature gradient between grain pairs that are unconnected by binder.In contrast to the decrease of heat flux on the faces for Neumann boundary conditions, Dirichlet boundary conditions induce a slight increase of heat flux on the faces, but only in places where the faces cut through a sand grain which has higher conductivity than its surroundings.The leftmost and rightmost cutouts in Figure 17 highlights a solitary sand grain just barely cutting the face of the unit cell.For Dirichlet boundary conditions, the temperature on the faces is fixed while the temperature gradient is applied on all interior points of the domain.For Neumann boundary conditions, the temperature field in this single grain is not limited and is therefore not subject to high heat fluxes.The fluctuation of the temperature field is shown in Figure 18.The different boundary conditions impact the fields of temperature fluctuations most, especially the zero values on Dirichlet boundary.Because the prescribed mean of the temperature gradient is the same as for the other boundary conditions, the field of temperature fluctuations shows larger domains of positive fluctuation on the right side of the unit cell and an inverse effect on the other side.For Neumann boundary conditions, the values of high temperature gradients on the bottom edge is reflected in the temperature fluctuation as well, where artifacts on the boundary occur, that we do not observe for periodic boundary conditions.
When comparing the computed effective properties, Schneider et al. 86 showed that the results do not perfectly agree with experimental results.An explanation for this disagreement considered by the authors is that the simulated boundary conditions do not match the experiment.Table 4 shows the computed effective conductivities for two samples of the sand grain microstructure with all three types of boundary conditions.We observe that, like for the stochastic microstructures considered before, Dirichlet boundary conditions lead to mismatch of the apparent properties obtained by experiment and simulation for both types of sand.Even though periodic boundary conditions lead to lower apparent values, the conductivities still exceed the experimental data.
With Neumann boundary conditions, we obtain apparent boundary conditions noticeably lower than periodic and Dirichlet.These apparent values are close to the experimental values.The Neumann boundary conditions mirror the reduced conductivity between the grains that Schneider et al. suspected to be the root cause for the reduced effective values during the experiment.The Dirichlet boundary conditions might keep the temperature fluctuation at the boundary stable, but cannot prevent the unchecked increase of temperature in the interior sand grains, leading to instability in the solver.Even though for infinitely large volume elements the apparent properties generally align with the effective properties, this realistic example shows that boundary conditions still influence obtained apparent properties significantly.Therefore, to approximate the effective properties by computing the apparent properties using a computational homogenization approach, it of utmost importance to know the actual boundary conditions with which the experimental results were produced.Then, enabled by the methods introduced in this work, the use of the suitable boundary conditions is possible, leading to a fast and efficient approximation of the microstructured material's apparent properties.

CONCLUSION
The work at hand was concerned with the efficient computational resolution of the cell problems arising in stochastic homogenization of thermal conductivity problems.More precisely, we considered discretizations on a regular grid and established a unified computational framework for Green's operator methods for Dirichlet, Neumann, and periodic boundary conditions.In particular, the unified perspective permitted to transfer results well-known for periodic boundary conditions, for example, bounds on the iteration count depending only on the initial condition of the iteration and the material contrast, to the other boundary conditions.Thus, Green's operator methods remove the dependence of the iteration count on the resolution which may be responsible for a weak performance of computational methods.
To evaluate the corresponding Green's operators, we followed two essentially different strategies.The first strategy is based on Fourier-related transforms, and exploits analytical expressions for the inverse of the finite-difference Laplacian in terms of the coefficients of the transform.Whereas this approach is standard for periodic boundary conditions via the discrete Fourier transform (DFT), the approach is less commonly employed for Dirichlet as well as Neumann boundary conditions utilizing the discrete sine and cosine transforms (DST/DCT), at least in the setting of heterogeneous coefficients.
The second strategy to obtain Green's operator is based on low-rank techniques.More precisely, although the finite-difference Laplacian may be represented as a sparse matrix, its inverse requires a fully populated matrix.However, this inverse may be approximated by a matrix with low ranks to high precision.For the work at hand, we made use of the TT format which turned out to be quite powerful when inverting the Laplacian in high dimensions. 43As we are mostly interested in two and three space dimensions, the QTT format turns out to be particularly favorable for low-rank techniques. 55We used a method based on exponential sums and the scaling and squaring algorithm to compute a low-rank QTT approximation of the Green's operators.Moreover, it might be interesting to investigate a third strategy, where the action of Green's operator is performed in real space. 42e implemented the outlined approaches for the EJ finite-volume discretization of thermal conductivity.Computational experiments clearly showed that, despite identical computational complexity, the practical performance of the Fourier-related transforms is significantly better than the low-rank strategies, outperforming the latter by at least an order of magnitude.The reasons for these differences are found, on the one hand, in the high implementation quality of the Fourier-related transforms, for example, by exploiting cache alignment and parallelization capabilities, and, on the other hand, in the fact that Green's operator is represented exactly by the Fourier-related transforms, whereas the use of a mere approximation to the Green's operator is detrimental to the performance of the QTT approach.Although (Q)TT strategies have their merits, in particular in high dimensions, they appear to be suboptimal for the application considered in this article, essentially due to the strong competition presented.
With these findings, we investigated the influence on the imposed boundary conditions on the computed apparent properties of materials with random microstructures.We studied three different ensembles of matrix-inclusion composites: perfectly periodic cells, extracted microstructures and those with hard boundaries, for example, where the inclusions are not permitted to intersect the boundary.It turned out that periodic boundary conditions (for the temperature-fluctuation field) worked nicely under all circumstances.Dirichlet boundary conditions led to a strong overestimation of the effective conductivities, and appeared to work reasonably well for hard boundaries only.In turn, this insight provides a motivation for engineers with a Dirichlet-only computational software at hand to exploit microstructure-generation tools with hard boundaries.Neumann boundary conditions were inferior to periodic boundary conditions, but provided a reasonable use case for engineering purposes.
Last but not least, we studied the influence of the boundary condition on the computed apparent properties of sand-grain microstructures.In fact, in case representative volume elements are excessively large, experimental tests provide only apparent properties, and the influence of the boundary conditions may be significant.Therefore, accounting for these conditions is absolutely critical for validating microstructure models.
The established methodology is expected to provide a common working ground for computational homogenization methods on regular grids.The framework is useful for other physical problems, as well, for example, problems in solid mechanics, and providing appropriate Green's operators for these scenarios is desirable.

2
The explicit jump grid for Dirichlet and Neumann boundary conditions for a two-dimensional domain Y 2 = [0, L] 2 with N = 4 voxels and an element size h = L∕N.The positions of temperature values are marked in red.The temperature gradient values on the voxel faces are marked in blue and green for the derivatives in x-and y-directions, respectively.Boundary values set to zero in round brackets.(A) Dirichlet boundary conditions.(B) Neumann boundary conditions.

F
I G U R E 11 Iteration count for a material ratio  = 10 2 and for radii 2r 1 = r 2 = 0.4.The FFT-based preconditioner on an N 3 EJ grid was used.F I G U R E 12Extracted microstructures with increasing edge length L = 16, 32, 64, … , 512 m and spheres of diameter 10 m.

F
I G U R E 13 Two-dimensional sketch of artificial microstructures with spherical inclusions showing different methods for generating the unit cell highlighted in gray.The microstructure is characterized by either (A) periodicity of the gray cell in all directions, (B) an extraction in x-direction with periodicity of the gray cell in the other directions, or (C) a hard wall in x-direction marked in blue with periodicity of the gray cell in the other directions.F I G U R E 14 Effective conductivity for microstructures from spherical inclusions with diameter 10 m with  = 1.5 W∕(K m) and a matrix with  = 0.195 W∕(K m).Standard deviation is denoted by the error bars.

F I G U R E 15
Sand grain microstructure B with edge length 1546 m and voxel size 3 m.The binder is highlighted in red.TA B L E 3Conductivity values in W/(K m) of the sand core components.85

F
I G U R E 16 x-component of the temperature gradient in K∕m on a slice of a sand grain microstructure with different boundary conditions for a uniaxial load in x-direction.(A) Neumann.(B) Periodic.(C) Dirichlet.F I G U R E 17 Heat flux in W∕m 2 on sand grain microstructure A with different boundary conditions.Uniaxial mean temperature gradient in x-direction.(A) Neumann.(B) Periodic.(C) Dirichlet.In Figure 17, the heat flux field for different boundary conditions is shown.Starting with the Neumann boundary condition, we observe some heterogeneous artifacts on the domain edges as the heat flux is forced to be zero on the domain edges.As the heat flux is decreased by the boundary condition near the edges of the domain, the solution field for periodic boundary conditions shows a zone of high heat flux between two grain cores that Neumann boundary conditions do not replicate.Like the temperature gradient fields, the boundary conditions do not influence the solution fields on the interior of the unit cell in any meaningful way.A slight decrease of the high flux zones where binder connects two sand grains can be observed for Neumann boundary conditions when compared to Dirichlet boundary conditions.

F I G U R E 18
Temperature fluctuation in K on sand grain microstructure A with different boundary conditions.(A) Neumann.(B) periodic.(C) Dirichlet.TA B L E 4 Computed effective conductivity in W/(K m) and experimental data from Reference 85 for microstructure samples A and B.