# Development and comparison of spectral algorithms for numerical modeling of the quasi-static mechanical behavior of inhomogeneous materials

In the current work, a number of algorithms are developed and compared for the numerical solution of periodic (quasi-static) linear elastic mechanical boundary-value problems (BVPs) based on two different discretizations of Fourier series. The first is standard and based on the trapezoidal approximation of the Fourier mode integral, resulting in trapezoidal discretization (TD) of the truncated Fourier series. Less standard is the second discretization based on piecewise-constant approximation of the Fourier mode integrand, yielding a piecewise-constant discretization (PCD) of this series. Employing these, fixed-point algorithms are formulated via Green-function preconditioning (GFP) and finite-difference discretization (of differential operators; FDD). In particular, in the context of PCD, this includes an algorithm based on the so-called "discrete Green operator" (DGO) recently introduced by Eloh et al. (2019), which employs GFP, but not FDD. For computational comparisons, the (classic) benchmark case of a cubic inclusion embedded in a matrix (e.g., Suquet, 1997; Willot, 2015) is employed. Both discontinuous and smooth transitions in elastic stiffness at the matrix-inclusion (MI) interface are considered. In the context of both TD and PCD, a number of GFP- and FDD-based algorithms are developed. Among these, one based on so-called averaged-forward-backward-differencing (AFB) is shown to result in the greatest improvement in convergence rate. As it turns out, AFB is equivalent to the "rotated scheme" (R) of Willot (2015) in the context of TD. In the context of PCD, comparison of the performance and convergence behavior of AFB/R- and DGO-based algorithms shows that the former is more efficient than the latter for larger phase contrasts.

There are no comments yet.

## Authors

• 1 publication
• 1 publication
• 1 publication
• 1 publication
âˆ™ 04/21/2020

### Gradient discretization of two-phase flows coupled with mechanical deformation in fractured porous media

We consider a two-phase Darcy flow in a fractured porous medium consisti...
âˆ™ 02/12/2021

### A space-time discretization of a nonlinear peridynamic model on a 2D lamina

Peridynamics is a nonlocal theory for dynamic fracture analysis consisti...
âˆ™ 09/20/2021

### Energy-stable discretization of two-phase flows in deformable porous media with frictional contact at matrix-fracture interfaces

We address the discretization of two-phase Darcy flows in a fractured an...
âˆ™ 11/11/2020

### Gradient discretization of two-phase poro-mechanical models with discontinuous pressures at matrix fracture interfaces

We consider a two-phase Darcy flow in a fractured and deformable porous ...
âˆ™ 05/16/2021

### Positivity preservation of implicit discretizations of the advection equation

We analyze, from the viewpoint of positivity preservation, certain discr...
âˆ™ 07/08/2021

### Monolithic multigrid for a reduced-quadrature discretization of poroelasticity

Advanced finite-element discretizations and preconditioners for models o...
âˆ™ 05/27/2019

### DBFFT: A displacement based FFT approach for non-linear homogenization of the mechanical behavior

Most of the FFT methods available for homogenization of the mechanical r...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 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

 f=Â¯f+~f,Â¯f:=1lâˆ«l0f(x)Â dx,~f:=fâˆ’Â¯f, (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

 u=uh+up,uh(x):=c+Â¯Hx,up(x):=âˆ«~H(x)Â dx, (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,

 ~E(x)Â¯E=~S(x)Â¯S,up(x)Â¯E=âˆ«~S(x)Â¯SÂ dx, (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

 Ï•Ïµ(x;c)=12+12tanh((xâˆ’c)/Ïµ) (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 (),

 Î½Ïµ(x):=Ï•(x;cl,Ïµ)âˆ’Ï•(x;cr,Ïµ) (5)

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

 S=(1âˆ’Î½Ïµ)SM+Î½ÏµSI=SM+SIMÎ½Ïµ,SIM:=SIâˆ’SM, (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

 Fmf(x):=âˆ‘mÎº=âˆ’meÄ±kÎºx^fÎº,^fÎº:=1lâˆ«l0eâˆ’Ä±kÎºxf(x)Â dx, (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

 ^ftÎº:=nâˆ’1/2Ë‡ftÎº,Ë‡ftÎº:=nâˆ’1/2âˆ‘ni=1eâˆ’Ä±kÎºxifti, (9)

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

 fT(x):=Ï†tÎ¼(x)Ë‡ftÎ¼,Ï†tÎ¼(x):=nâˆ’1/2eqÎ¼x,fti=Ft+Î¼iË‡ftÎ¼,Ë‡ftÎ¼=Ftâˆ’Î¼ifti, (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,

 upT(x)=Ï†Î¼(x)Ë‡upÎ¼,âˆ‡aupT(x)=Ï†Î¼(x)Ë‡upÎ¼âŠ—qaÎ¼, (11)

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

 âˆ‡FDupT(x):=hâˆ’1[upT(x+h)âˆ’upT(x)]=Ï†Î¼(x)Ë‡upÎ¼âŠ—qFDÎ¼=:~EFD(x),âˆ‡BDupT(x):=hâˆ’1[upT(x)âˆ’upT(xâˆ’h)]=Ï†Î¼(x)Ë‡upÎ¼âŠ—qBDÎ¼=:~EBD(x),âˆ‡CDupT(x):=12hâˆ’1[upT(x+h)âˆ’upT(xâˆ’h)]=Ï†Î¼(x)Ë‡upÎ¼âŠ—qCDÎ¼=:~ECD(x),âˆ‡hCupT(x):=hâˆ’1[upT(x+12h)âˆ’upT(xâˆ’12h))]=Ï†Î¼(x)Ë‡upÎ¼âŠ—qhCÎ¼=:~EhC(x), (12)

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

 qFDÎ¼:=hâˆ’1(eqÎ¼hâˆ’1),qBDÎ¼:=hâˆ’1(1âˆ’eâˆ’qÎ¼h), (13)

and

 qCDÎ¼:=hâˆ’1sinh(qÎ¼h),qhCÎ¼:=2hâˆ’1sinh(12qÎ¼h). (14)

Analogous to (11), TD of yields

 TT(x)=Ï†Î¼(x)Â Ë‡TÎ¼,divbTT(x)=Ï†Î¼(x)Â Ë‡TÎ¼qbÎ¼, (15)

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

 âˆ’divbCH[âˆ‡aupT(x)]=Ï†Î¼(x)Â Ë‡Gâˆ’1H(qaÎ¼,qbÎ¼)Ë‡upÎ¼ (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.,

 qbÎ¼=qaâˆ—Î¼ (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

 âˆ‡BDf(x+12h)=âˆ‡hCf(x) (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)

 fPC(x):=Ï†cÎ¼(x)Ë‡fcÎ¼,Ï†cÎ¼(x):=scÎ¼Ï†tÎ¼(x),fci=Fc+Î¼iË‡fcÎ¼,Ë‡fcÎ¼=Fcâˆ’Î¼ifci, (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

 upPC(x)=Ï†cÎ¼(x)Â Ë‡upcÎ¼,TPC(x)=Ï†cÎ¼(x)Â Ë‡TcÎ¼, (22)

and

 âˆ‡aupPC(x)=Ï†cÎ¼(x)Â Ë‡upcÎ¼âŠ—qaÎ¼,Â Â divbTcPC(x)=Ï†cÎ¼(x)Ë‡TcÎ¼qbÎ¼. (23)

Further,

 âˆ’divbCH[âˆ‡aupPC(x)]=Ï†cÎ¼(x)Â Ë‡Gâˆ’1H(qaÎ¼,qbÎ¼)Ë‡upcÎ¼, (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.

As evident, this algorithm is independent of from (20). This is in contrast to the PCD-based algorithm of Eloh etÂ al. (2019), to which we now turn.

### 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

 Fcpf(x):=âˆ‘pâˆ’1Îº=âˆ’peÄ±kÎºx^fcÎº=nâˆ’1/2âˆ‘nâˆ’1Ï‰=0eÄ±kÏ‰xâˆ‘mâˆ’1Î½=âˆ’meÄ±kÎ½nxsÎ½n+Ï‰Ë‡fcÎ½n+Ï‰ (25)

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

 Ë‡fcÎ½n+Ï‰=(âˆ’1)Î½Ë‡fcÏ‰,nâˆ’1/2âˆ‘nâˆ’1Ï‰=0eÄ±kÏ‰xciË‡fcÏ‰=f(xci)=fci, (26)

from (20) results in the PCD

 fEl(x):=Ï†ElÏ‰(x)Ë‡fcÏ‰,Ï†ElÏ‰(x):=nâˆ’1/2âˆ‘mâˆ’1Î½=âˆ’meÄ±kÎ½n+Ï‰x(âˆ’1)Î½sÎ½n+Ï‰, (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 transform222To 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.

 ^Î“H(k)^T(k)=âˆ’^GH(k)^T(k)kâŠ—k=0,kâ‰ 0, (28)

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

Here, and