## 1 Introduction

Numerical modeling of the mechanical behaviour of inhomogeneous materials is often based on the assumption that their microstructure is periodic. In this case, a well-known alternative to finite difference, finite element, or isogeometric, methods for the numerical solution of the corresponding boundary value problems (BVPs) is offered by spectral methods (e.g., Gottlieb and Orszag, 1977; Canuto etÂ al., 1988; Trefethen, 2000; Press etÂ al., 2007; Kopriva, 2009). Besides the possibility of exponential convergence, spectral methods generally offer much higher spatial resolution than finite-difference- or finite-element-based ones. Indeed, as discussed for example by Press etÂ al. (2007, Â§20.7), when applicable, these are the methods of choice for this purpose. Another difference here is that, in contrast to finite-difference and finite-element methods, which approximate / discretize the model field relations, spectral methods approximate / discretize their solution. In particular, in the Fourier case, this latter is based on discretization of truncated Fourier series, referred to as Fourier series discretization in what follows.

For the current case of numerical solution of the (quasi-static) momentum balance and corresponding mechanical BVPs, Fourier series discretization has traditionally been combined with differential-operator preconditioning (e.g., Moulinec and Suquet, 1994; Suquet, 1997)

, formally analogous to differential-operator inversion for analytic solution of BVPs in mathematical physics based on Green functions (e.g., Lippmann-Schwinger). This is referred to as "Green-function preconditioning" (GFP) in what follows. Together with continuous and discrete Fourier transformation, such methods have a long history of application in material science and in particular in continuum micromechanics

(e.g., Khachaturyan, 1983; Mura, 1987; Moulinec and Suquet, 1994; Suquet, 1997; Chen, 2002) to the modeling of material microstructure. The corresponding mechanical BVP is formulated at the level of a unit cell or representative volume element. Traditionally, these are strain-based, geometrically linear, and based on fixed-point iterative solution (e.g., Suquet, 1997; Eyre and Milton, 1999; Michel etÂ al., 2000, 2001; Lebensohn, 2001; Eisenlohr etÂ al., 2013; Lebensohn and Rollett, 2020). More recently, these have been extended conjugate-gradient- and Newton-Krylov-based numerical solution and geometric nonlinearity (e.g., Kabel etÂ al., 2014; Shanthraj etÂ al., 2015; Mishra etÂ al., 2016).As well-known, in the case of material heterogeneity due to discontinuous material properties, Fourier series discretizationn suffers from oscillations due to the Gibbs phenomena and to aliasing (e.g., Trefethen, 2000; Kopriva, 2009). In general, both have an adverse effect on convergence behavior, resulting in a significant loss of algorithmic efficiency and robustness. To address these issues, a number of algorithms going beyond the original "basic scheme" of Suquet (1997) have been developed. One class of such algorithms combines Fourier series discretization with finite-difference discretization (FDD) of differential operators (e.g., Dreyer and MÃ¼ller, 2000; Brown etÂ al., 2002; Moulinec and Silva, 2014; Willot, 2015), resulting in a reduction of the effect of oscillations on convergence behavior due to modes at wavelengths shorter than the nodal or grid spacing. As shown in particular by Willot (2015), among such "accelerated schemes" (e.g., Moulinec and Silva, 2014), his "rotated scheme" is most effective in improving algorithmic efficiency and robustness. This scheme has been employed for numerical solution of BVPs for heterogeneous materials in a number of works (e.g., for field dislocation mechanics in Djaka etÂ al., 2017). More recently, Schneider etÂ al. (2017) showed that finite-element discretization based on linear hexahedral elements with reduced integration is equivalent to the rotated scheme of Willot (2015). Quite recently, Lucarini and Segurado (2019) developed a displacement-based approach via direct Fourier transformation and real-function-based Fourier transform symmetry reduction, yielding an algorithmic system based on a full-rank associated Hermitian matrix amenable to preconditioning. For the latter purpose, they worked with a Green function based on the average elastic stiffness. In contrast to Willot (2015) and the current work, FDD was not employed.

All algorithms or "schemes" discussed up to this point employ Fourier series discretization in "standard" form, i.e., based on trapezoidal approximation / discretization (TD) of the integral for Fourier modes with respect to the unit cell. TD is the basis of Fourier interpolation and the standard form of the discrete Fourier transform

(e.g., Trefethen, 2000; Press etÂ al., 2007; Kopriva, 2009). Alternatively, Brisard and Dormieux (2010) and Anglin etÂ al. (2014) discretize by assuming the solution is piecewise-constant, resulting in piecewise-constant discretization (PCD). This approach has been exploited recently by Eloh etÂ al. (2019), who combine PCD with GFP to obtain an algorithm depending on a so-called discrete Green operator (DGO) for numerical solution of periodic mechanical BVPs for materially heterogeneous elastostatics via fixed-point interation.One purpose of the current work is to develop additional algorithms in the context of Fourier series discretization based on TD and PCD. All of these exploit GFP, and some of them FDD. A second purpose is to carry out detailed theoretical and computational comparisons of these with existing algorithms, in particular with the "rotated" (R) scheme of Willot (2015), and with the DGO-based scheme of Eloh etÂ al. (2019). The corresponding computational comparisons are based on the classic benchmark case of a periodic matrix-inclusion microstructure with cubic inclusions having sharp corners. The material inhomogeneity takes the form of (both discontinuous and smooth) phase contrast in elastic stiffness (e.g., Suquet, 1997; Nemat-Nasser and Hori, 1999; Willot, 2015) across the matrix-inclusion (MI) interface. In contrast, Eloh etÂ al. (2019) focus mainly on material heterogeneity due to eigenstrains (see also e.g., Anglin etÂ al., 2014); in one example, however, spherical inclusions are considered. Likewise, Lucarini and Segurado (2019) focus (solely) on spherical inclusions.

To establish relevant basic concepts and issues, the work begins in Section 2 with basic considerations in one dimension (1D). These include analytic solution of the periodic matrix-inclusion mechanical BVP and its truncated Fourier series representation. Attention is focused next on the development of algorithms in 1D based on TD and PCD of truncated Fourier series in Section 3. As discussed above, all of these employ GFP, and some of them FDD. Computational comparison of these is then carried out in Section 4

for the matrix-inclusion benchmark case. In particular, material heterogeneity in the form of both (i) discontinuous and (ii) smooth compliance / stiffness distributions are considered. Multidimensional forms of the 1D algorithms are obtained via direct tensor product generalization in Section

5. Corresponding computational comparsions are presented in Section 6. The work ends with a summary and discussion in Section 7.In this work, (three-dimensional) Euclidean vectors are represented by lower-case bold italic characters

. In particular, let represent the Cartesian basis vectors. Second-order tensors are represented by upper-case bold italic characters , with the second-order identity. By definition, any maps any linearly into a vector . Fourth-order Euclidean tensors are denoted by upper-case slanted sans-serif characters. By definition, any maps any linearly into a second-order tensor . Let (summation convention) represent the scalar product of two arbitrary-order tensors and . Using this, defines the transpose of . Let represent the symmetric part of . The tensor products and will also be employed. Additional notation will be introduced as needed in what follows.## 2 Basic considerations in 1D

### 2.1 Analytic solution of periodic boundary-value problem

Let the interval of length represent the unit cell of the periodic matrix-inclusion (MI) microstructure. In what follows, we work with the split

(1) |

of any (integrable) into mean and "fluctuating" parts.

Let represent the displacement field on the unit cell and the corresponding distortion. Integration of yields the split

(2) |

of into "homogeneous" and "particular" parts, with the constant of integration. In what follows, and are assumed given or known. Then is determined, and represents the primary unknown field.

Excluding cracking, kinematic compatibility requires , and so , to be continuous everywhere, in particular at MI interfaces. In addition, quasi-static mechanical equilibrium requires the stress to be continuous everywhere and in fact constant in 1D; then and . In this case, the linear elastic relation for the 1D strain in terms of the compliance and stiffness results in the equilibrium relations and for in 1D. In turn,

(3) |

then hold for and , the latter via (2).

In the classic composite case (e.g., Suquet, 1997), the interface between two bulk phases is idealized as "sharp" in the sense that material properties like are assumed to vary discontinuously across the interface. On the other hand, phase field models idealize such interfaces as a mixture of the bulk phases in which properties like vary smoothly from one bulk phase to the other. Assuming the classic case of a double well potential and "flat" interface region of half-width centered at , solution of the equilibrium Ginzburg-Landau relation at zero stress yields

(4) |

for the phase field of the MI system varying between (matrix) and (inclusion) in the MI interface region. The "sharp interface" limit of this is given by in terms of the modified Heaviside step function (i.e., , , ; e.g., Bracewell, 2000, Chapter 4). Assuming that the left (right) MI interface is at (),

(5) |

represents the 1D "volume density" of the inclusion. In terms of this,

(6) |

holds for the elastic compliance of the MI unit cell (note ). Substituting (6) into (3) yields

(7) |

with and . Note that has been assumed in (7) without loss of physical generality. and are displayed in Figure 1.

As expected from (3), determines a continuous and smooth for , and a discontinuous for . Likewise, is continuous and smooth for , but only continuous for .

### 2.2 Approximation based on truncated Fourier series

As a first step toward Fourier-series-based numerical solution of the above BVP, consider next the truncated Fourier series approximation of and from (7). The truncated Fourier series of any is given by

(8) |

(e.g., Trefethen, 2000; Kopriva, 2009; Liu, 2011) with and . As is well-known, in general at one or more . In particular, for discontinuous on , then, deviates from due to (i) at points of discontinuity , and (ii) truncation of . Since (i) is related to the Gibbs phenomenon, it will be referred to as Gibbs error in what follows. In the sharp () MI interface case, then, is affected by (i) and (ii), whereas is affected only by (ii). These are displayed in Figure 2.

As expected, at ; rather, . Being due to approximation of discontinuous by , note that Gibbs error is unaffected by truncation of to , i.e., independent of . In addition, it is independent of the magnitude of .

## 3 Numerical solution algorithms in 1D

### 3.1 Algorithms based on trapezoidal discretization

As discussed in the introduction, algorithms for the numerical solution of the current BVP are based on (i) discretization of (8) in 1D, and (ii) tensor-product generalization of these to 3D. More specifically, two discretizations of (8) are considered and compared in this work. The most common of these is trapezoidal discretization of (8) based on that of with nodes at and nodal spacing . In this case, (8) and discretizes to

(9) |

("t" stands for "trapezoidal").
As detailed in Appendix A, (9)
induces the trapezoidal discretization
(TD)^{1}^{1}1In (10) and
all relations to follows, repeated indices, e.g.,
or , imply a sum over these from to
when this makes sense.

(10) |

of (8) via (A.6) for . Here, , , , and from (A.7). If is prescribed, note that is determined for a given discretization. For notational simplicity, the "" superscript for "trapezoidal" is dropped in what follows.

As discussed in the introduction, following previous work (e.g., Dreyer and MÃ¼ller, 2000; Moulinec and Silva, 2014; Willot, 2015), finite difference discretization (FDD) of differential operators in the context of (10) is employed here for improved robustness in convergence. For example,

(11) |

via FDD of in terms of the effective wavenumber . In particular,

(12) |

for forward (), backward (), central (), and half central (), difference discretization, respectively, with

(13) |

and

(14) |

Analogous to (11), TD of yields

(15) |

via FDD of the divergence operator. Together, (11) and (15) induce the preconditioned form

(16) |

of the homogenous elliptic operator with

(17) |

the corresponding Green function ().

In what follows, is determined from a given choice of in two different ways. The first is based on conjugacy, i.e.,

(18) |

(e.g., Willot, 2015, Equations (19) and (24)). For example, for from (13), or for via (14). The second way is based on the choice and the observation that determines via the mean-value theorem. The choice and identity

(19) |

then result in the combination . For this choice, then, both the stress divergence and displacement are calculated at . Note that this choice does not satisfy (18).

The above results are incorporated in Algorithm TD1.

The direct Fourier case is included here by defining ("F" for "Fourier"). Rather than being strain-based like the "basic scheme" of Suquet (1997) and Michel etÂ al. (2000, 2001), note that Algorithm TD1 is displacement-based. As evident from (17), is not invertible when either or vanish. For even, (13) and (14) imply at for , and at for . For this latter case, we follow Willot (2015, Equation (29)) and formally set for . For odd, only at for all cases.

### 3.2 Algorithms based on piecewise-constant discretization

Assume now that in (8) is constant on each subinterval of with value at the mid-point or "center" of . In this case, in (8) reduces to

(20) |

with the (unnormalized) sinc function. As shown by comparison of with from (9), differences between these include (i) the dependence of on the sinc "filter" , and (ii) different discretizations of . Leaving the details to Appendix B, (20) results in the piecewise-constant discretization (PCD)

(21) |

of (8) via (B.8) for , analogous to (10) in the TD case, with and . In particular, and . Analogous to (11) and (15), then, we have

(22) |

and

(23) |

Further,

(24) |

analogous to (16). Completely analogous to the case of Algorithm TD1, then, these relations and the fact that can be employed to obtain Algorithm PCD1.

### 3.3 Algorithm of Eloh etÂ al. (2019)

Following Brisard and Dormieux (2010), Eloh etÂ al. (2019) recently developed an algorithm also based on PCD and (20) quite different than Algorithm PCD1. To this end, they formulate their algorithm based on from (8) for and (20), and truncate the result. Equivalently, as done here, one can work with the truncated form

(25) |

of from the start (for even). Combining this with (the second of)

(26) |

from (20) results in the PCD

(27) |

(sum over repeated ) of (8)
for ,
analogous to (21) ("El" stands for "Eloh").
As in the basic scheme (e.g., Suquet, 1997), they work with (i) strain
as the primary discretant, and
(ii) the Fourier transform^{2}^{2}2To be more precise, again like
Suquet (1997), Eloh etÂ al. (2019) work with the
polarization stress instead of
directly as done in Algorithm DGO1. This is also true for
the 3D case and Algorithm DGO2 below.

(28) |

of mechanical equilbrium in pre-conditioned Lippmann-Schwinger form. Again leaving the details to Appendix B, this results in Algorithm DGO1.

Here, and

Comments

There are no comments yet.