1 Introduction↩︎

1.1 State of the art↩︎

Accurate comparisons between experimental results and simulations requires robust and efficient numerical methods capable of handling large 3D grids. Since the seminal work of [1], [2], fast Fourier transform (FFT)-based solvers establish themselves as efficient alternatives to standard finite element (FE) codes, in describing the macroscopic properties of heterogeneous media, as well as the distribution of local fields (e.g., strain, stress). For an extensive overview of classical FFT-based solvers and their applications, we refer the reader to the review articles [3]–[5].

FFT-based solvers exploit the application of a discrete Green operator in Fourier space to solve the governing equations efficiently. The initial proposition of [2] can be regarded as a collocation method, starting from the continuous equation and using a truncation of Fourier series to solve the so-called Lippman-Schwinger equation for an auxiliary problem defined on a homogeneous material. Another possibility consists of using Galerkin variational formulations [6]–[8] to build the discrete Green operator. Moreover, using a finite difference scheme to obtain the discrete differential operators (e.g., divergence, gradient), a discrete Green operator can directly be built from the local equations defined on discrete grids. In the last three decades, the improvements in the expression of the discrete Green operator coupled with acceleration techniques [3], [9], [10] of the so-called original basic scheme [1], [2], have allowed to study more complex diffusion and mechanical problems. In this context, the present work proposes a new variant of the Anderson acceleration based on the divergence of the polarization stress field.

Despite their effectiveness and growing popularity, classical FFT-based solvers remain inherently limited to periodic media and periodic boundary conditions (BCs), which restricts their applicability to a broader range of problems. In the recent years, various propositions have been discussed to circumvent this limitation. A buffer zone (with an arbitrary behavior) in which the desired heterogeneous medium is embedded, is used to mimic Neumann and/or Dirichlet BCs [9], [11], [12]. This approach has been successfully applied for linear and non-linear material behaviors [13], [14], at a minimal cost in classical FFT-based solvers. Nevertheless, the buffer zone is shown to influence the convergence rate [13]. Furthermore, studying an effective conductivity problem, [15] employed the mirror transformation method to apply uniform BCs. With this approach, the problem is solved on computation domain 8 times larger than the original desired 3D unit-cell, leading to an increase in the simulation time. [10] proposed a displacement based boundary domain integral equation that enables the solution of elasticity boundary value problems. The applicability of this approach is demonstrated on 2D elastic composites subjected to homogeneous loading.

Another strategy to handle non-periodic BCs with FFT-based solvers consists of using discrete trigonometric transforms [16], [17]. With this method the desired problem (i.e., partial differential equations) is solved without enlarging the domain. It relies on employing sine and cosine real transforms instead of the discrete Fourier transform (DFT) as in the classical FFT-based solvers. In this sense, [17] develop limited number of mixed uniform BCs for isotropic linear elastic materials. [18], [19] also apply this approach to represent only Dirichlet BCs in a deformation-based formulation. Using a Galerkin-based method to build the discrete Green operator, [20] were able to discuss heterogeneous elastic structures subjected to non-periodic mixed Dirichlet/Neumann BCs. Meanwhile, [21], [22] treated static and transient diffusion problems build on finite difference scheme approach in deriving of the discrete Green operator. More recently, damage issues in 2D composites are discussed in the work of [23] under non-periodic BCs (Dirichlet and free surfaces). The DTTs-based approach for applying non-periodic BCs, is the subject of the present work. We aim to provide a robust and generic FFT-based solver for accounting any possible BCs (non-periodic Neumann-Dirichlet and/or a mix with Periodic), implemented within a parallel framework (distributed memory) and suitable for small and finite transformations.

As stated earlier, the discrete Green operator can be evaluated with various approaches [3], [8]. In this work, we build on the finite difference methods. In the context of periodic condition, [24] proposed a rotated scheme (hereafter referred to as hexahedral scheme in light of FE [25]) helping to achieve in general case, better convergence in comparison with the basic scheme of Moulinec and Suquet [1], [2]. In general, within the fixed-point method, the hexahedral scheme still suffers convergence issue with microstructures that include void phases [26]. Very recently, keeping with periodic BCs, [27] proposed a double tetrahedron finite difference scheme that remedied these convergence issues. With the latter scheme, stress and strain fields exhibiting fewer numerical artifacts can also be obtained. Staggered grid scheme used in the work of [26] provide ringing-free local fields and good convergence in the context of Dirichlet BCs. The effort is made in this work to compare the double tetrahedron scheme to the hexahedral scheme in the context of non-periodic BCs, and particularly in non-linear conditions.

1.2 Contributions↩︎

In this context, we propose an efficient, versatile and robust parallel FFT-based solver that seamlessly accommodates arbitrary combinations of periodic and/or non-periodic BCs in small and finite transformation simulations of the linear and non-linear mechanical responses of small-scale heterogeneous structures, relying on a formulation in displacement.

For the purpose, the mechanical problem for Cauchy continuum theories is formulated in a displacement-based approach. By introducing the fluctuation displacement as an unknown, the general (non-periodic and/or periodic) BCs are linked with the symmetries and/or periodicity of the relevant fluctuation displacement and stress components. We take advantage of the fact that for a given field, symmetry conditions allow to compute the DTTs which themselves are analytically related to the DFT of the “virtually” extended periodic field. The application of the discrete Green operator in the novel generic FFT-based solver can then be interpreted as an appropriate use of the classical periodic discrete Green operator on the later “virtually” extended periodic field. Particularities introduced by non-periodic BCs in the expressions of the discrete Green operators are also discussed based on the hexahedral [24] and double tetrahedron [27] finite difference schemes. We show that the classical periodic discrete Green operator based on a homogeneous isotropic elastic medium restricts the total number of applicable BCs. Reliable alternatives [19], [20] are discussed. Sec. 2 details the major steps in building the generic FFT-based solver.

The implementation is performed in a parallel framework in the open-source AMITEX\(^\star\) solver (new version, in progress, of the AMITEX_FTTP code [28]). The parallel implementation relies on a 2D pencil decomposition of the unit-cell that allows for distributing the simulation over a larger number of processors than a simple 1D slab decomposition, originally available in the FFTW library [29] (a 1D slab decomposition of a \((n^v)^3\) unit-cell is limited to \(n^v\) processors, with \(n^v\) the number of voxel in a given direction). This 2D decomposition is managed in AMITEX\(^\star\) by the open-source library 2DECOMP&FFT [30], [31]. Moreover, DTTs and DFT (collectively referred to as discrete transforms, DTs) are computed using the FFTW library [29] through the interface 2DECOMP&FFT. A key contribution of the present work is the extension of 2DECOMP&FFT to support DTTs and a mix with DFT, which were not available in the original library [30], [31].

The purpose of Sec. 3 is to demonstrate the robustness and versatility of the solver using non-trivial problems, with various BCs that could not be taken into account by classical FFT-based solvers, linear and non-linear complex behaviors, small and finite transformations (including void growth and coalescence). The advantage of the 2D decomposition is also demonstrated in a crystal plasticity simulation. Finally, concluding remarks and recommendations are presented in Sec. 4.

Notation↩︎

Unless otherwise indicated, the Einstein index convention is used. Volume mean values are denoted \(\left< \star \right>\); discrete Fourier transform: \(\widehat{\star} = \mathrm{DFT}(\star)\); discrete transforms: \(\widetilde{\star} = \mathrm{DT}(\star)\) accounting more generally for trigonometric and Fourier transforms; vector \(\boldsymbol{a}\); second-order tensor \(\oalign{\boldsymbol{A}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\); deviatoric part of a second-order tensor \(\mathrm{dev}(\oalign{\boldsymbol{A}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth})\); fourth-order tensor \(\mathbb{A}\); dyadic product “\(\otimes\)” such as \(\left(\boldsymbol{a} \otimes \boldsymbol{b}\right)_{ij} = a_ib_j\); Hadamard product “\(\odot\)” such as \(\left(\oalign{\boldsymbol{A}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\odot\oalign{\boldsymbol{B}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\right)_{ij} = A_{ij}B_{ij}\) (no summation on the indices). \((\boldsymbol{e}_1,\boldsymbol{e}_2,\boldsymbol{e}_3)\) represents the classical Cartesian basis of the euclidean space and \((O,\boldsymbol{e}_1,\boldsymbol{e}_2,\boldsymbol{e}_3)\) denote the orthonormal coordinate system.

2 Generic FFT-based parallel solver for non-periodic and/or periodic BCs↩︎

As aforementioned, the classical FFT-based solver is inherently restricted to the scope of periodic BCs. The objective of this section is to motivate and detail the different modifications associated with non-periodic BCs. To proceed, we hereby present the studied domain and its discretization on a regular grid. This enables us to rewrite the system of equations to be solved for the mechanical problem (based on Cauchy continuum theories). An iterative (accelerated) displacement-based fixed-point method is used. We then describe the trigonometric transforms used with the aim to preserve the structure of the accelerated FFT-based solver. Additionally, the finite difference options employed to build the discrete Green operators are outlined. This section is concluded with implementation specifications related to parallel computing and the treatment of combined periodic and non-periodic BCs.

2.1 Domain discretization↩︎

In the context of FFT-based solver, the domain of interest on which the mechanical problem is to be solved, is considered as a parallelepiped domain \(\Omega\) as shown in Fig. 1 (a). The domain \(\Omega\) is delimited by faces denoted by \(S_{i\alpha}\) with \(i\) being the direction of the normal to the given face and \(\alpha \in \{0,1\}\) the position of the face (\(\alpha=0\) corresponds to the face that contains the origin point at the coordinates \(\left[0,0,0\right]\)). \(\boldsymbol{n}_{i\alpha}\) is the outer normal to the face \(S_{i\alpha}\).

The continuum domain \(\Omega\) is discretized (i.e., \(\Omega^d\)) by elementary voxels, of sizes \(h_i,\,i \in \{1,2,3\}\) (in each \(\boldsymbol{e}_i\) direction, in this work, we use \(h_1=h_2=h_3\)). For a given direction \(\boldsymbol{e}_i\), the number of considered voxels is labeled \(n_i^v\) (leading to a total of \(n_1^v \times n_2^v \times n_3^v\) voxels or \(N_1^n\times N_2^n \times N_3^n\) nodes with \(N_i^n=n_i^v+1\)). Each voxel is characterized by its eight nodes with Cartesian coordinates \(\left[j_1 h_1, j_2 h_2, j_3 h_3\right]\) with \(j_i\in \llbracket 0,n_i^v\rrbracket,\,i \in \{1,2,3\}\). Subsequently, the voxel centers are positioned at the coordinates \(\left[(j_1+ 1/2) h_1, (j_2+ 1/2) h_2, (j_3+ 1/2) h_3\right]\) with \(j_i\in \llbracket 0,n_i^v-1\rrbracket,\,i \in \{1,2,3\}\).

The degrees of freedom of the mechanical problem (typically the displacement vector field) are classicaly defined on the nodes. According to the finite difference schemes used (which are presented later), the strains and stresses (second-order tensorial fields), which are inputs/outputs of a given behavior law \(\mathcal{F}\), are defined at voxel centers. Fig. 1 (b) illustrates an example of a \(4\times5\) voxels centers (or alternatively \(5\times6\) nodes) discretization. This is a 2D projection of the voxels connected to the face \(S_{10}\) (i.e., voxels centers and nodes are not located on the same plane).

2.2 Formulation of the mechanical problem↩︎

To define the mechanical problem, we introduce the following notations and physical fields. \(\mathcal{F}\) is the user-defined behavior law, \(\boldsymbol{u}\) denotes the total displacement vector field, \(\boldsymbol{u^*}\) is an arbitrary displacement vector field which values at the faces \(S_{i\alpha}\) correspond to the desired applied Dirichlet BCs. The fluctuation displacement field [18], [20] is then defined as \(\boldsymbol{u^f} = \boldsymbol{u} - \boldsymbol{u^*}\). An applied stress traction vector field \(\boldsymbol{t^*}\) (or \(\boldsymbol{T^*}\) in finite transformation framework) which values at the faces \(S_{i\alpha}\) corresponds to the desired applied Neumann BCs, \(\oalign{\boldsymbol{\sigma}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\) is the Cauchy second-order stress tensor, \(\oalign{\boldsymbol{P}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\) is the first Piola-Kirchoff second-order stress tensor, \(\oalign{\boldsymbol{\Pi}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\) corresponds to the second Piola-Kirchoff second-order stress tensor. \(\oalign{\boldsymbol{P}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\) and \(\oalign{\boldsymbol{\Pi}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\) are related to \(\oalign{\boldsymbol{\sigma}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\) through the gradient of the transformation \(\oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \oalign{\boldsymbol{1}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} + \oalign{\boldsymbol{\nabla u}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\), according to the formulas \[\oalign{\boldsymbol{P}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \mathrm{det}(\oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}) \, \oalign{\boldsymbol{\sigma}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} \oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^{-T} \quad;\quad\oalign{\boldsymbol{\Pi}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \mathrm{det}(\oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}) \, \oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^{-1} \oalign{\boldsymbol{\sigma}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} \oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^{-T}\]

In the context of general BCs, it is possible to impose a given BC type (periodic, Dirichlet or Neumann) per face \(S_{i\alpha},\,i \in \{1,2,3\}, \alpha \in \{0,1\}\) and per component of the displacement or traction vector field. Moreover, small and finite transformation frameworks can be distinguished when formulating the boundary value problem in mechanics. In fact, with small strain hypothesis, the system of equations is given by \[\text{Unknown} \,\,\boldsymbol{u} / \left\{ \begin{array}{l} \boldsymbol{u} = \boldsymbol{u^f} + \boldsymbol{u^*}\\[0.5em] \boldsymbol{\nabla}\cdot\oalign{\boldsymbol{\sigma}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} =0 \quad \text{in} \quad \Omega \\[0.5em] \oalign{\boldsymbol{\sigma}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \mathcal{F} ( \oalign{\boldsymbol{\nabla^s u}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}) \quad \text{in} \quad \Omega \\[0.5em] \text{On each face } S_{i\alpha},\,i \in \{1,2,3\}, \alpha \in \{0,1\} \text{ and for each component } q \in \{1,2,3\} \\[0.5em] \qquad \bullet \quad \text{ if } \text{ Periodic BC: } u^f_{q} \text{ periodic and } \sigma_{iq} \text{ periodic} \\[0.5em] \qquad \bullet \text{ or if } \text{ Dirichlet BC: } u_{q} = u_{q}^{*} \, \\[0.5em] \qquad \bullet \text{ or if } \text{ Neumann BC: } \sigma_{iq} \left(-1\right)^{\alpha+1}=t^*_{q} \end{array} \right.\] We recall that, for a given field, periodic condition means that the field has the same value on opposite faces and anti-periodic condition expresses the fact that the considered field has opposite values on opposite faces. For this reason, the stress component \(\sigma_{iq}\) being periodic is equivalent to term \(\left[\oalign{\boldsymbol{\sigma}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\,\boldsymbol{n}_{i\alpha}\right]_{q}\) being anti-periodic. Within finite transformation framework, we use a Lagrangian formulation so that the system of equation is now \[\text{Unknown} \,\,\boldsymbol{u} / \left\{ \begin{array}{l} \boldsymbol{u} = \boldsymbol{u^f} + \boldsymbol{u^*}\\[0.5em] \boldsymbol{\nabla}\cdot\oalign{\boldsymbol{P}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} =0 \quad \text{in} \quad \Omega \\[0.5em] \oalign{\boldsymbol{P}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \mathcal{F} ( \oalign{\boldsymbol{\nabla u}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}) \quad \text{in} \quad \Omega \\[0.5em] \text{On each face } S_{i\alpha},\,i \in \{1,2,3\}, \alpha \in \{0,1\} \text{ and for each component } q \in \{1,2,3\} \\[0.5em] \qquad \bullet \quad \text{ if } \text{ Periodic BC: } u^f_{q} \text{ periodic and } P_{iq} \text{ periodic} \\[0.5em] \qquad \bullet \text{ or if } \text{ Dirichlet BC: } u_{q} = u_{q}^{*} \, \\[0.5em] \qquad \bullet \text{ or if } \text{ Neumann BC: } P_{iq} \left(-1\right)^{\alpha + 1}=T^*_{q} \end{array} \right. \label{eq:TF95initial}\tag{1}\] As the systems of equations in both small and finite transformation cases are similar, we then focus on the most general not approximated formulation of finite transformation in the following analysis. These continuous systems are solved on the discretized domain \(\Omega^d\), replacing continuous derivation operators by finite difference derivation operators. For the purpose, we recall that displacement fields are defined on the voxel nodes and strain/stress fields on voxel centers.

Observing the well-known particular case of periodic BCs, the periodic extension of displacement and stress fields around the concerned faces, allows to compute the gradient of the displacement and the divergence of the stress in the vicinity of the faces. Actually, to evaluate the displacement gradient at voxel centers connected to the unit-cell faces, using a given finite difference scheme, the displacements at nodes located on the faces are defined by periodicity. Correspondingly, the computation of the stress divergence at nodes located on the faces, using a given finite difference scheme, requires the stress values at mirror voxel centers located outside the computational domain, which are obtained by periodic extension. Following the same idea, non-periodic BCs will rely on the definition of appropriate extensions. To this end, the symmetry (even-symmetry) and anti-symmetry (odd-symmetry) are used. A given field that possesses an anti-symmetric extension w.r.t. to a face, has a null value on the concerned face. In contrast, a symmetric extension does not prescribe any specific value on the face, leaving the corresponding quantity unconstrained.

From the definition of the fluctuation displacement field, \(\boldsymbol{u^{f}} = \boldsymbol{u} - \boldsymbol{u^*}\), in the case of Dirichlet BC of the component \(q \in \{1,2,3\}\), on a given face \(S_{i\alpha},\,i \in \{1,2,3\}, \alpha \in \{0,1\}\), the condition \(u_{q} = u_{q}^{*}\) then yields to a vanishing fluctuation displacement field of the corresponding component, i.e., \(u_q^{f}=0\) on \(S_{i\alpha}\). Therefore, Dirichlet BC can be associated with an anti-symmetric extension of the fluctuation field component \(u_q^{f}\) w.r.t. to the face \(S_{i\alpha}\). In addition, in order to always satisfy the corresponding equilibrium condition (i.e., \(\left[{}_d\boldsymbol{\nabla}\cdot\oalign{\boldsymbol{P}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\right]_q =0\)) at the face \(S_{i\alpha}\), symmetric extension of the component \(P_{iq}\) and anti-symmetric extension of the other stress components \(P_{il}, l \neq q\), are used. Note that for the component \(P_{iq}\), the symmetric extension let the quantity free (unconstrained) at the corresponding face.

In case of Neumann BC of the component \(q \in \{1,2,3\}\), on a given face \(S_{i\alpha},\,i \in \{1,2,3\}, \alpha \in \{0,1\}\), as the stress fields are classically defined on the voxel centers, the Neumann BC, \(P_{iq}\left(-1\right)^{\alpha+1}=T^*_{q}\), is not automatically enforced at the face \(S_{i\alpha}\) (voxels faces). In order to illustrate the analysis on this case and without any loss of generality, we consider an 1D configuration as described in Fig. 2. The desired stress field \(P\) in the domain \(\Omega^d\) (represented in solid black line in Fig. 2) possesses a Neumann BC, i.e., \(P=T^*\), only on the right face (for simplicity purpose). Similarly to the Dirichlet BCs, an applied stress field \(P^*\) is introduced such as \(P^*\) is null at all voxel centers inside the discretized domain \(\Omega^d\) and non-null with a value \(P^*=2T^*\) at the “virtual” mirror voxel center w.r.t. to the faces \(S_{i\alpha}\) (blue line in Fig. 2). With this choice, the value of \(P^*\) interpolated at the face \(S_{i\alpha}\) is equal to \(T^*\). A fluctuation stress field \(P^f\) can be defined as \(P^f = P-P^*\). By construction, the fluctuation stress field strictly coincides with the desired field \(P\) inside the unit-cell \(\Omega\) and vanishes at the face \(S_{i\alpha}\). This is in accordance with the 3D case for the hexahedral and double tetrahedron finite difference schemes that will be discussed in this work. One can conclude that a Neumann BC of a component \(q \in \{1,2,3\}\) w.r.t. to a face \(S_{i\alpha},\,i \in \{1,2,3\}, \alpha \in \{0,1\}\) is characterized by an anti-symmetric extension of the corresponding fluctuation stress component \(P^f_{iq}\). The remaining components \(P^f_{il}, l\neq q\) are considered to have no particular fixed values (even-symmetry extension w.r.t. to the face \(S_{i\alpha}\)). The displacement component \(u_q=u^f_q\) to be determined as consequence of Neumann BC is then also assumed to be symmetric w.r.t. to the face \(S_{i\alpha}\).

Given the established connection between the non-periodic BCs and the symmetries of the fields, the discrete form (depending on the discretization and a finite difference form) of the boundary value problem in Eq. 1 , is expressed as follows \[\text{Unknown} \,\,\boldsymbol{u^f} / \left\{ \begin{array}{l} \boldsymbol{u} = \boldsymbol{u^f} + \boldsymbol{u^*}\\[0.5em] \oalign{\boldsymbol{P}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \oalign{\boldsymbol{P^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} + \oalign{\boldsymbol{P^*}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} \\[0.5em] {}_d\boldsymbol{\nabla}\cdot\oalign{\boldsymbol{P}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} =0 \quad \text{in} \quad \Omega^d \\[0.5em] \oalign{\boldsymbol{P}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \mathcal{F} ( {}_d\oalign{\boldsymbol{\nabla u}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}) \quad \text{in} \quad \Omega^d \\[0.5em] \text{On each face } S_{i\alpha},\,i \in \{1,2,3\}, \alpha \in \{0,1\} \text{ and for each component } q \in \{1,2,3\} \\[0.5em] \quad \bullet \quad \text{ if } \text{ Periodic BC: } u^f_{q} \text{ periodic and } P_{iq} \text{ periodic} \\[1em] \quad \bullet \text{ or if } \text{ Dirichlet BC: } u_{q}^{f} \text{ anti-symmetric, } P_{iq} \text{ symmetric and } P_{il} \text{ with } l\neq q \text{ anti-symmetric } \\[2em] \quad \bullet \text{ or if } \text{ Neumann BC: } u_{q}^{f} \text{ symmetric, } P^f_{iq} \text{ anti-symmetric and } P^f_{il} \text{ with } l\neq q \text{ symmetric } \end{array} \right. \label{eq:GeneralBCtobesolved}\tag{2}\] In practice, the components of \(\oalign{\boldsymbol{P}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\) (or, equivalently, \(\oalign{\boldsymbol{P}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^f\)) exhibit the same symmetries and/or periodicity as the components of \(\oalign{\boldsymbol{\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\). This prompts an implementation by specifying only the spatial extensions of the fluctuation displacement vector \(\boldsymbol{u^f}\) for all general BCs (periodic, Dirichlet or Neumann). Subsequently, based on the extensions of the fluctuation displacement fields in the cases of Dirichlet and Neumann BCs, the standard FFT-based iterative solver [2], [24] can not be directly applied. Adaptations are required for a generic FFT-based solver that can be used for all general BCs.

2.3 Overview of the displacement-based iterative fixed-point approach↩︎

In order to clarify the novel features of the FFT-based solver in the context of the non-periodic BCs, we first review the general displacement-based iterative fixed-point approach used in this work.

Knowing that FFT-based solvers are classically used in the context of periodic BCs, the mechanical problem to solve is reduced to the following one \[\left\{ \begin{array}{l} \boldsymbol{u} = \boldsymbol{u^f} + \boldsymbol{u^*}\\[0.5em] {}_d\boldsymbol{\nabla}\cdot\oalign{\boldsymbol{P}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} =0 \quad \text{in} \quad \Omega^d \\[0.5em] \oalign{\boldsymbol{P}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \mathcal{F} ( {}_d\oalign{\boldsymbol{\nabla u}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}) \quad \text{in} \quad \Omega^d \\[0.5em] \text{On all faces } S_{i\alpha},\,i \in \{1,2,3\}, \alpha \in \{0,1\} \text{ and for all components } q \in \{1,2,3\} \\[0.5em] \qquad \bullet \text{ Periodic BC: } u^f_{q} \text{ periodic and } P_{iq} \text{ periodic} \end{array} \right. \label{eq:fullperiodic}\tag{3}\] A priori pre-conditioner \(\mathbb{R}_0:\oalign{\boldsymbol{\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\) is introduced so that the mechanical problem in Eq. 3 is rewritten as \[\left\{ \begin{array}{l} \boldsymbol{u} = \boldsymbol{u^f} + \boldsymbol{u^*}\\[0.5em] {}_d\boldsymbol{\nabla}\cdot \left(\mathbb{R}_0:{}_d\oalign{\boldsymbol{\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\right) = {}_d\boldsymbol{\nabla}\cdot(\mathbb{R}_0:{}_d\oalign{\boldsymbol{\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} - \oalign{\boldsymbol{P}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}) \mathrel{\vcenter{:}}= \boldsymbol{p} \quad \text{in} \quad \Omega^d \\[0.5em] \oalign{\boldsymbol{P}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \mathcal{F} ( {}_d\oalign{\boldsymbol{\nabla u}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}) \quad \text{in} \quad \Omega^d \\[0.5em] \text{On all faces } S_{i\alpha},\,i \in \{1,2,3\}, \alpha \in \{0,1\} \text{ and for all components } q \in \{1,2,3\} \\[0.5em] \qquad \bullet \text{ Periodic BC: } u^f_{q} \text{ periodic and } P_{iq} \text{ periodic} \end{array} \right.\] \(\oalign{\boldsymbol{\tau}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \mathbb{R}_0:{}_d\oalign{\boldsymbol{\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} - \oalign{\boldsymbol{P}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\) is the so-called polarization stress [2]. The principle of the fixed-point algorithm is to predict \(\left[\boldsymbol{u^f}\right]^{m+1}\) the fluctuation displacement field at an iteration \(m+1\) given \(\left[\boldsymbol{u^f}\right]^{m}\) at the iteration \(m\). First, we express, at the iteration \(m\), the divergence of the polarization stress \(\boldsymbol{p}^{m}\) \[\boldsymbol{p}^{m} = {}_d\boldsymbol{\nabla}\cdot\left(\mathbb{R}_0:\left[{}_d\oalign{\boldsymbol{\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\right]^{m} - \left[\mathcal{F} ( {}_d\oalign{\boldsymbol{\nabla u}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth})\right]^{m}\right) \label{eq:div95pol}\tag{4}\] which is used in a second step to deduce \(\left[\boldsymbol{u^f}\right]^{m+1}\) for the iteration \(m+1\), based on the expression \[{}_d\boldsymbol{\nabla}\cdot \left(\mathbb{R}_0:\left[{}_d\oalign{\boldsymbol{\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\right]^{m+1}\right) = \boldsymbol{p}^{m} \label{eq:div95pol95next95iter}\tag{5}\] The relationship in Eq. 4 is evaluated in real space and the last operation in Eq. 5 can be evaluated in Fourier space using DFT, resulting in the following expression \[\widehat{\left[\boldsymbol{u^f}\right]^{m+1}} = {}_d\widehat{\mathrm{DGO}}(\widehat{\boldsymbol{p}^{m}}) \label{eq:GreenOp}\tag{6}\] defining the Fourier discrete Green operator \({}_d\widehat{\mathrm{DGO}}\) (core of the FFT-based solver), which depends a priori on the finite difference scheme and where \(\widehat{G}[k] = \displaystyle \sum_{j=0}^{n_{points}-1} g[j] \, \exp\left(-\mathrm{i} \,\frac{2 \pi k}{n_{points}} j\right)\) for a field \(G\) defined on \(n_{points}\) points (for simplicity, a 1D example is used knowing that a 3D transform is a succession of 1D transform). Finally, \(\left[\boldsymbol{u^f}\right]^{m+1}\) is obtained by an inverse DFT calculation. \(\left[\boldsymbol{u^f}\right]^{m+1}\) is then reintroduced in Eq. 4 for a new iteration of the fixed-point algorithm up to the equilibrium convergence criterion, expressed as follows (for cuboid voxels, see [13] for general non-cubic voxels) \[\sqrt{\frac{\|{}_d\boldsymbol{\nabla}\cdot\oalign{\boldsymbol{P}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\|_{L^2}}{\|\oalign{\boldsymbol{P}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\|_{L^2} } } \,h_1 < \epsilon \label{eq:equilibrium95condition}\tag{7}\] with \(\epsilon\) the tolerance. It is well established that this fixed-point approach is efficient in solving the problem in Eq. 3 . In order to preserve this efficiency in solving the general BCs as described in Eq. 2 , it necessitates to be able to compute DFT and the discrete Green operator in general cases (possibility of mixing Dirichlet, Neumann and/or periodic BCs). The subsequent sections address these points.

In addition, it is important to note that the basic iterative algorithm can be accelerated by different approaches [3], [5]. In this work, an Anderson convergence acceleration technique is used [32], [33]. The idea is directly taken from the acceleration technique used in the FE code CAST3M [34] to accelerate a modified Newton-Raphson algorithm. Indeed, every 2 or 3 iterations, a displacement field \(\left[\boldsymbol{u^f}\right]^{m+1}\) for the iteration \(m+1\) is proposed based on the obtained fields, \(\left[\boldsymbol{u^f}\right]^{j}\), and residual fields, \(\left[\boldsymbol{u^f}\right]^{j} - \left[\boldsymbol{u^f}\right]^{j-1}\), at the previous \(4\) couples \(j\)-iterations. This technique described in [9], has been shown to significantly enhance the convergence rate of the basic scheme in classical FFT-based solver.

2.4 Generic FFT-based solver features: linking DTTs and DFT↩︎

In this section, we build on the link established in the previous Sec. 2.2 between the general BCs and the symmetry and/or periodicity extensions of the desired fields to develop the generic FFT-based solver. This approach aims to facilitate the preservation of the classical FFT-based iterative fixed-point algorithm, with certain modifications, in order to address the problem formulated in Eq. 2 (possibility to combine periodic BCs with non-periodic BCs). Depending on the BCs on each face, symmetry (even-symmetry), anti-symmetry (odd-symmetry) or periodicity extensions are proposed for each component of the fluctuation displacement and stress fields. According to these extensions, the fluctuation displacement and stress fields can be “virtually” extended to obtain periodic fields from which the discrete periodic Green operator can be build (see Eq. 6 ). This approach requires to express the relation between the DTs of the original fields and the DFT of the “virtually” extended fields. For the sake of clarity, without any loss of generality, we can restrict to a 1D scalar field defined at the nodes (we recall that to the displacement field \(u^f\) is defined at the nodes).

We consider a discretized 1D scalar field \(G\) with values \(g[j],\,j\in\llbracket 0,n^v \rrbracket\) at the nodes (\(N^n=n^v+1\) nodes in total) which possesses periodicity or symmetry extensions on respectively the left (\(j=0\)) and right (\(j=n^v\)) sides. In the case of periodic condition, the “virtually” extended field \(E\) is then considered to be strictly equal to the primary field \(G\). In the non-periodic case, four possibilities can be distinguished for the node-defined field. For a field \(G\) with symmetry extensions, the original field can be “virtually” extended either twice or four times (symmetry factor \(S_f=2\) or \(S_f=4\)) in order to obtain a “virtual” extended periodic field \(E\). Fig. 3 illustrates this concept through two scenarios: a field \(G\) with anti-symmetric/anti-symmetric extensions (i.e., Dirichlet/Dirichlet BCs) on both sides and another field \(G\) with anti-symmetric/symmetric extensions (i.e., Dirichlet/Neumann BCs) on the left and right sides. On the extended periodic field \(E\), a computation of a DFT, \(\widehat{E}\), can be performed. It can be demonstrated that there is a connection between the DFT of \(E\), \(\mathrm{DFT}(E)[k]=\widehat{E}[k] = \displaystyle \sum_{j=0}^{S_f(N^n-1)-1}e[j] \exp\left(-\mathrm{i} \dfrac{2\pi k}{S_f(N^n-1)} j\right)\), and the DTs of \(G\), \(\mathrm{DT}(G)=\widetilde{G}\). Calculation of the DFT and the DTs is provided in Tab. 1 and specially, their links are summarized in the third column of Tab. 1 depending on the BCs. It is shown that for a field, with a given BC, once the DT of the original field is computed on the desired unit-cell (see index \(j\) in the summation sign in Tab. 1 restricted to the domain of interest), the DFT of the extended field is obtained using the constant value (real or purely imaginary), denoted \(\alpha^{\mathrm{BC}}\) in Tab. 1. Furthermore, it is worth noting that the value of fundamental frequencies, denoted by \(\xi^k\) with \(k\) being the index of selected frequencies and which are used in the calculations of the DTs, \(\widetilde{G}\), correspond to appropriately selected frequencies for the DFT, \(\widehat{E}\) (see the third column of Tab. 1). The non-null stored values in Fourier space depend on the BCs. For the generic AMITEX\(^\star\) solver, DTT calculations are performed in the FFTW library [29] (which also manage the DFT operations), relying on new upgrades of the open-source library 2DECOMP&FFT [30], [31] as an interface.

To sum up, even in the presence of non-periodic BCs, the FFT-based solver (mainly the application of the discrete Green operator, see Eq. 6 ) can be used by capitalizing on the relationships between DTs and DFT. In this sense, in order to harmonize the treatment of all possible BCs, we decide to use a \(4\) times extended “virtual” domain for fields that possess any type of symmetry or periodicity extensions (i.e., \(S_f=4\) for all cases).

Alg. 4 summarizes the different steps of the resolution with the generic FFT-based solver. Moreover, it can be seen in Alg. 4 that Anderson acceleration is applied to \(\boldsymbol{p}\), the divergence of the polarization field. To the best of the authors’ knowledge, this constitutes a novel departure from the conventional approach [9], in which Anderson acceleration is applied directly to the fluctuation displacement field, \(\boldsymbol{u^f}\). Our numerical simulations indicate that the \(\boldsymbol{p}\)-based formulation often yields better convergence performance than the conventional \(\boldsymbol{u^f}\)-based approach.

Going forward, it is still crucial to elucidate the symmetry and/or periodicity extensions of the divergence of the polarization stress, \(\boldsymbol{p}\), as well as the methodology for deriving the discrete Green operator based on the finite difference scheme and the BCs. To do so, the derivation operations present in the model, here the gradient, the divergence and the Laplacian (combination of the divergence and the gradient), are first expressed in the next section.

2.5 Finite difference schemes↩︎

In the mechanical problem expressed in Eq. 2 , discrete gradient operations are imperative to quantify strains at the voxel centers and divergence operations to test the equilibrium equation at the nodes. The present work uses finite difference schemes to evaluate the discrete derivation operations. Due to their demonstrated advantages in the context of periodic BCs, the hexahedral scheme [24], [25] and the recently proposed double tetrahedron scheme [22], [27], are investigated in the present paper. The flexibility of the present implementation in the AMITEX\(^\star\) solver allows for other types of finite difference schemes or finite element types.

To ensure the completeness of the present work, the fundamental concepts of the schemes are reviewed here, prior to the analysis of the modifications introduced by non-periodic BCs. Fig. 5 illustrates the two finite difference schemes. The hexahedral scheme (see Fig. 5 (a)) can be viewed as a linear hexahedral FE with a single “integration point” [25], and the double tetrahedron scheme (see Fig. 5 (b)) can also be regarded, from an implementation point of view, as a linear hexahedral FE but with two “integration points” [35]. This indicates that, for both schemes, the derivation of a node-defined field (e.g., a displacement field) results in a field defined at the voxel center (e.g., a strain field), and vice versa. The hexahedral (HEX1) scheme uses all eight points (voxel nodes or centers) altogether in a single derivation operation. Meanwhile, for the double tetrahedron (TETRA2) case, if the field is initially defined at the voxel nodes, then two derivatives are obtained, each of them being computed with one of the two regular tetrahedron defined in the voxel (four nodes are collectively involved at once). When the field is defined at centers (with two values per voxel, one per “integration point”), the resulting derivative is defined at nodes (with a single value per node).

The accelerated fixed-point FFT-based solver, necessitates the computation of the gradient, divergence operators in real and Fourier spaces to deduce the discrete Green operator. For the purpose, the following two sections detail the derivation operations in real and Fourier spaces for both finite difference schemes. A cross validation of the derivatives in real and Fourier spaces is performed, for any combination of BCs in the Sec. 2.7.3. Unless otherwise stated and without any loss of generality, we considered a 3D scalar field \(G\) defined at the voxel nodes (identified by \(j_i \in \llbracket 0, n_i^v \rrbracket,\,i \in \{1,2,3\}\)) and with general BCs in a given direction w.r.t. to a face. As stated before, the first derivative is obtained at the voxel centers (voxel-defined field).

2.5.1 Hexahedral scheme↩︎

In this section, the HEX1 finite difference scheme [24] is introduced in real and Fourier spaces (depending on the BCs).

• Derivation in real space:

For the 3D node-defined scalar field \(G\), we can write down the component of the derivative in each direction at the voxel centers (identified by \(j_i \text{ with } j_i + 1/2 \in \llbracket 0, n_i^v-1 \rrbracket,\,i \in \{1,2,3\}\)) \[\begin{align}D_1^HG\left[j_1+\frac{1}{2},j_2+\frac{1}{2},j_3+\frac{1}{2}\right] = \displaystyle \frac{1}{4h_1} \Bigl[ & G[j_1+1,j_2,j_3] - G[j_1,j_2,j_3] \nonumber \\ + & G[j_1+1,j_2+1,j_3] - G[j_1,j_2+1,j_3] \nonumber \\ + & G[j_1+1,j_2,j_3+1] - G[j_1,j_2,j_3+1] \nonumber \\ + & G[j_1+1,j_2+1,j_3+1] - G[j_1,j_2+1,j_3+1] \Bigr] \label{eq:D1HG} \end{align}\tag{8}\] \[\begin{align} D_2^HG\left[j_1+\frac{1}{2},j_2+\frac{1}{2},j_3+\frac{1}{2}\right] = \displaystyle \frac{1}{4h_2} \Bigl[ & G[j_1,j_2+1,j_3] - G[j_1,j_2,j_3] \nonumber \\ + & G[j_1+1,j_2+1,j_3] - G[j_1+1,j_2,j_3] \nonumber \\ + & G[j_1,j_2+1,j_3+1] - G[j_1,j_2,j_3+1] \nonumber \\ + & G[j_1+1,j_2+1,j_3+1] - G[j_1+1,j_2,j_3+1] \Bigr] \label{eq:D2HG} \end{align}\tag{9}\] \[\begin{align} D_3^HG\left[j_1+\frac{1}{2},j_2+\frac{1}{2},j_3+\frac{1}{2}\right] = \displaystyle \frac{1}{4h_3} \Bigl[ & G[j_1,j_2,j_3+1] - G[j_1,j_2,j_3] \nonumber \\ + & G[j_1+1,j_2,j_3+1] - G[j_1+1,j_2,j_3] \nonumber\\ + & G[j_1,j_2+1,j_3+1] - G[j_1,j_2+1,j_3] \nonumber \\ + & G[j_1+1,j_2+1,j_3+1] - G[j_1+1,j_2+1,j_3] \Bigr] \label{eq:D3HG} \end{align}\tag{10}\]

For the mechanical problem in this work, it is interesting to express the gradient and Laplacian (gradient followed by divergence) of a node-defined displacement vector field \(\boldsymbol{v}\) (with scalar components \(v_j,\, j \in \{1,2,3\}\)) and the divergence of a voxel-defined second-order stress/strain tensor \(\oalign{\boldsymbol{S}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\) (with scalar components \(S_{ij},\, i,j \in \{1,2,3\}\)). In this sense, we have \[\oalign{\boldsymbol{{}_{H}\nabla v}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = D_i^Hv_j \,\boldsymbol{e}_i \otimes \boldsymbol{e}_j \quad; \quad \boldsymbol{{}_{H}\nabla}\cdot \oalign{\boldsymbol{S}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = D_j^HS_{ij} \,\boldsymbol{e}_i \quad; \quad \boldsymbol{{}_{H}\Delta v} = \boldsymbol{{}_{H}\nabla} \cdot \oalign{\boldsymbol{{}_{H}\nabla v}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = D_i^H D_i^Hv_{j} \,\boldsymbol{e}_j\]

• Derivative symmetry and/or periodicity extensions:

In general, given the symmetry and/or periodicity extensions of a field w.r.t. a face, the extensions of its derivatives can be determined. We consider a field \(G\) that possesses symmetry and/or periodicity extensions, denoted by the triplet \(ABC\), where \(A\), \(B\) and \(C\) are respectively the spatial extension in the three directions \(\boldsymbol{e}_1\), \(\boldsymbol{e}_2\) and \(\boldsymbol{e}_3\) (i.e., w.r.t. the faces \(S_{1\alpha}\), \(S_{2\alpha}\) and \(S_{3\alpha}\), \(\alpha \in \{0,1\}\)). \(A\), \(B\) and \(C\) being a couple of conditions XY where X, Y are respectively the BCs at faces with \(\alpha=0\) (left side) and \(\alpha=1\) (right side). X, Y can be the periodicity (P), even-symmetry (S) or odd-symmetry (A): X and Y \(\in \left\{\text{P},\text{A},\text{S}\right\}\). If \(\text{X} \equiv \text{A}\), we note \(\overline{\text{X}}\equiv \text{S}\) and vice versa. If \(\text{X} \equiv \text{P}\) then \(\text{Y} \equiv \text{P}\) and \(\overline{\text{X}} \equiv \overline{\text{Y}} \equiv \text{P}\). Using the notation \(A^{'}B^{'}C^{'}\) for the extensions of the first-order continuous derivatives \(D_{i}G\), \(i \in \{1,2,3\}\), we have \[\left\{ \begin{array}{cccccc} \text{if } i=1: & A^{'}= \overline{A} &\text{and} & B^{'}= B &\text{and} & C^{'}= C \\[0.2em] \text{if } i=2: & A^{'}= A &\text{and} & B^{'}= \overline{B} &\text{and} & C^{'}= C \\[0.2em] \text{if } i=3: & A^{'}= A &\text{and} & B^{'}= B &\text{and} & C^{'}= \overline{C} \end{array} \right. \label{eq:Continuous95der95symmetry95per}\tag{11}\] For the HEX1-based discrete derivation, the parallelepiped voxel being itself invariant by symmetry and/or periodicity w.r.t. the faces \(S_{i\alpha}, \, i \in \{1,2,3\}\), the extensions of the discrete first-order derivative, \(D^H_{i}G\), at the faces, are consistent with the continuous analysis in Eq. 11 . This aspect is different for TETRA2 scheme.

• Derivation in Fourier space:

For any combination of BCs applied to a field \(G\), Eqs. 8 10 can be written in the Fourier space by applying DFT on the \(4\) times “virtually” extended field, \(E\). For simplicity of the notation and knowing the established link between the transforms \(\widehat{E}\) and \(\widetilde{G}\) (see Tab. 1 in Sec. 2.4), we use the notation \(\widehat{G}\) to signify the DFT of the “virtually” extended field based on the symmetries of \(G\). Therefore, combining Eqs. 8 10 with the DTT definitions introduced in Tab. 1, after elementary trigonometric manipulations, the DFT of the derivative components are expressed as \[\widehat{D_i^{H}\,G} =\mathrm{i} \, {}_{H}\xi_i \, \widehat{G} \, H_{tc} \qquad, \qquad i \in \{1,2,3\} \label{eq:derive95Fourier95HEX1}\tag{12}\] with \({}_{H}\xi_i , \, i \in \{1,2,3\}\) the components of the real value modified frequency vector, defined as \[{}_{H}\boldsymbol{\xi} = \begin{pmatrix} \displaystyle \frac{2}{h_1}\, \mathrm{sin}\left(\frac{\xi_1^{\mathrm{BCs}}}{2}\right)\,\mathrm{cos}\left(\frac{\xi_2^{\mathrm{BCs}}}{2}\right)\,\mathrm{cos}\left(\frac{\xi_3^{\mathrm{BCs}}}{2}\right) \\[1em] \displaystyle \frac{2}{h_2}\, \mathrm{cos}\left(\frac{\xi_1^{\mathrm{BCs}}}{2}\right)\,\mathrm{sin}\left(\frac{\xi_2^{\mathrm{BCs}}}{2}\right)\,\mathrm{cos}\left(\frac{\xi_3^{\mathrm{BCs}}}{2}\right) \\[1em] \displaystyle \frac{2}{h_3}\, \mathrm{cos}\left(\frac{\xi_1^{\mathrm{BCs}}}{2}\right)\,\mathrm{cos}\left(\frac{\xi_2^{\mathrm{BCs}}}{2}\right)\,\mathrm{sin}\left(\frac{\xi_3^{\mathrm{BCs}}}{2}\right) \end{pmatrix} \label{eq:modified95freq95hex1}\tag{13}\] and \[H_{tc} = \exp \left(\mathrm{i} \, \frac{\xi_1^{\mathrm{BCs}} + \xi_2^{\mathrm{BCs}} + \xi_3^{\mathrm{BCs}}}{2}\right) \label{eq:Htc95frequencies}\tag{14}\] The frequencies \(\xi_i^{\mathrm{BCs}},\, i \in \{1,2,3\}\) are the fundamental ones shown in Tab. 1 (i.e., variables \(\xi^k\) in Tab. 1), in accordance with the symmetry or periodicity extensions of the field \(G\). In the Fourier space, a derivative in one direction is then affected by the symmetry and/or periodicity extensions in the other directions. Moreover, when deriving in the Fourier space a field \(G\) defined at the voxel centers, the term \(H_{tc}\) in the Eq. 12 has to be replaced by its conjugate \(\overline{H_{tc}}\).

Similarly to the real space case, the gradient and Laplacian of a node-defined displacement vector field \(\boldsymbol{v}\) and, the divergence of a voxel-defined second-order stress/strain tensor \(\oalign{\boldsymbol{S}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\) are expressed in the Fourier space as \[\widehat{\oalign{\boldsymbol{{}_{H}\nabla v}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}} = \mathrm{i} \, \widehat{\boldsymbol{v}} \otimes {}_{H}\boldsymbol{\xi} \, H_{tc}\quad; \quad \widehat{\boldsymbol{{}_{H}\nabla}\cdot \oalign{\boldsymbol{S}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}} = \mathrm{i} \, \widehat{\oalign{\boldsymbol{S}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}}\,{}_{H}\boldsymbol{\xi} \, \overline{H_{tc}}\quad; \quad \widehat{\oalign{\boldsymbol{{}_{H}\Delta v}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}} = -\lVert {}_{H}\boldsymbol{\xi} \rVert^2 \widehat{\boldsymbol{v}} \label{eq:derivation95hex195fourier}\tag{15}\]

2.5.2 Double tetrahedron scheme↩︎

The finite difference scheme TETRA2 [27] is examined in the present section.

• Derivation in real space:

Defined at the voxel centers (identified by \(j_i + 1/2 \text{ with } j_i \in \llbracket 0, n_i^v-1 \rrbracket,\,i \in \{1,2,3\}\)), the derivatives in a chosen direction w.r.t. the regular tetrahedrons \(T_1\) and \(T_2\) (derivation supports), can be expressed as \[\left\{ \begin{array}{l} D_1^{T_1}G = \displaystyle \frac{1}{2h_1} \Bigl[ G[j_1+1,j_2+1,j_3] - G[j_1,j_2+1,j_3+1] + G[j_1+1,j_2,j_3+1] - G[j_1,j_2,j_3]\Bigr] \\[0.5cm] D_1^{T_2}G = \displaystyle \frac{1}{2h_1} \Bigl[ G[j_1+1,j_2,j_3] - G[j_1,j_2,j_3+1] + G[j_1+1,j_2+1,j_3+1] - G[j_1,j_2+1,j_3]\Bigr] \end{array} \right. \label{eq:D1T1T2G}\tag{16}\] \[\left\{ \begin{array}{l} D_2^{T_1}G = \displaystyle \frac{1}{2h_2} \Bigl[ G[j_1+1,j_2+1,j_3] - G[j_1,j_2,j_3] + G[j_1,j_2+1,j_3+1] - G[j_1+1,j_2,j_3+1]\Bigr] \\[0.5cm] D_2^{T_2}G = \displaystyle \frac{1}{2h_2} \Bigl[ G[j_1,j_2+1,j_3] - G[j_1+1,j_2,j_3] + G[j_1+1,j_2+1,j_3+1] - G[j_1,j_2,j_3+1]\Bigr] \end{array} \right. \label{eq:D2T1T2G}\tag{17}\] \[\left\{ \begin{array}{l} D_3^{T_1}G = \displaystyle \frac{1}{2h_3} \Bigl[ G[j_1,j_2+1,j_3+1] - G[j_1,j_2,j_3] + G[j_1+1,j_2,j_3+1] - G[j_1+1,j_2+1,j_3]\Bigr] \\[0.5cm] D_3^{T_2}G = \displaystyle \frac{1}{2h_3} \Bigl[ G[j_1,j_2,j_3+1] - G[j_1,j_2+1,j_3] + G[j_1+1,j_2+1,j_3+1] - G[j_1+1,j_2,j_3]\Bigr] \end{array} \right. \label{eq:D3T1T2G}\tag{18}\] Remark 1 is still valid. Additionally, for the finite difference scheme TETRA2, the second-order derivatives are calculated by combining both tetrahedrons (cross derivations), i.e., \(D_p^{T_1}\left(D_q^{T_2}G\right), \, p,q\in \{1,2,3\}\) or \(D_p^{T_2}\left(D_q^{T_1}G\right), \, p,q\in \{1,2,3\}\).

In an analogous manner to the HEX1 scheme, the gradient and the Laplacian of a node-defined displacement vector field \(\boldsymbol{v}\) (with scalar components \(v_j,\, j \in \{1,2,3\}\)) and the divergence of a voxel-defined second-order stress/strain tensor \(\oalign{\boldsymbol{S}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\) (with scalar components \(S_{ij},\, i,j \in \{1,2,3\}\)) are given in the form \[\oalign{\boldsymbol{{}_{T_r}\nabla v}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = D_i^{T_r}v_j \,\boldsymbol{e}_i \otimes \boldsymbol{e}_j\quad; \quad \boldsymbol{{}_{T_r}\nabla}\cdot \oalign{\boldsymbol{S}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = D_j^{T_r}S_{ij} \,\boldsymbol{e}_i \quad \text{with } r \in \{1,2\}\] \[\boldsymbol{{}_{T_1T_2}\Delta v} = \displaystyle \frac{1}{2} \left( D_i^{T_1} D_i^{T_2}v_{j} \, + D_i^{T_2} D_i^{T_1} v_{j} \,\right) \boldsymbol{e}_j\]

• Derivative symmetry and/or periodicity extensions:

In the previous case of the HEX1 scheme, the extensions of a derivative field are shown to be consistent with the continuous analysis. However, for the TETRA2 scheme, the symmetry and/or periodicity analysis are less straightforward. By construction, the tetrahedron \(T_2\) (respectively \(T_1\)) is the symmetry of the tetrahedron \(T_1\) (respectively \(T_2\)) w.r.t. the voxel faces (see Fig. 5 (b)). This means that the first-order derivation operation on a given tetrahedron is not invariant by symmetry. Considering a single tetrahedron, the first-order derivative of a field extended with appropriate symmetries do not recover any specific symmetries. In order to clarify this assertion, we consider the face \(S_{10}\) normal to the direction \(\boldsymbol{e}_1\) and containing here the node \(\left[0,0,0\right]\) (see Fig. 5 (b) with \(\left[j_1,j_2,j_3\right] = \left[0,0,0\right]\)). For the purpose of this demonstration, \(G\) is a field that possesses an even-symmetry w.r.t. the face \(S_{10}\). Combining the definition of the derivative in Eq. 16 \(_1\) and the symmetry condition of \(G\) w.r.t. the face \(S_{10}\), it follows that the derivative \(D_1^{T_1}G\) at the mirror voxel center is given by the relationship \[D_1^{T_1}G\left[j_1-\frac{1}{2},j_2+\frac{1}{2},j_3+\frac{1}{2}\right] = - D_1^{T_2}G\left[j_1+\frac{1}{2},j_2+\frac{1}{2},j_3+\frac{1}{2}\right]\] This expression shows no obvious symmetry property of the first-order derivative field \(D_1^{T_1}G\) w.r.t. the face \(S_{10}\) (the same holds true for the derivative field \(D_1^{T_2}G\)). The symmetries appear when considering a change of derivation support from \(T_1\) to \(T_2\). Meanwhile, for periodic BCs, both first-order derivative fields \(D_i^{T_1}G\) and \(D_i^{T_2}G\), \(i\in \{1,2,3\}\), are periodic. It is important to note that for the second-order derivatives, i.e., \(D_p^{T_1}\left(D_q^{T_2}G\right), \, p,q\in \{1,2,3\}\) or \(D_p^{T_2}\left(D_q^{T_1}G\right), \, p,q\in \{1,2,3\}\), the derived discrete fields have the same spatial support and also the same symmetry and/or periodicity as the field \(G\).

This discussion has a major consequence in the computation of the DFT of the first-order derivatives, in Fourier space. Actually, for the non-periodic case, the fact that the first-order derivative fields \(D_i^{T_r}G\), \(i \in \{1,2,3\}, r \in \{1,2\}\) have no specific symmetries, makes it impossible to use DTTs (on the original field to be derived) to evaluate these first-order derivatives in Fourier space. Nevertheless, the DTTs can be employed in the evaluation of the second-order derivatives, which are essential for the discrete Green operators. These issues do not emerge for full periodic BCs.

• Derivation in Fourier space:

Given the absence of discernible symmetry extensions of the first-order derivatives, it is not possible to evaluate these derivatives using DTTs. However, in the periodic case, since differentiation preserves periodicity, the DFT remains applicable. For the purpose, considering a periodic field, we have the DFT of the first-order derivatives as \[\begin{align} &\widehat{D_m^{T_1}\,G} = \mathrm{i} \, {}_{T_1}\xi_m \, \widehat{G} \, H_{tc} \quad; \quad m \in \{1,2,3\} \tag{19}\\ &\widehat{D_m^{T_2}\,G} = \mathrm{i} \, {}_{T_2}\xi_m \, \widehat{G} \, H_{tc} \quad; \quad m \in \{1,2,3\} \tag{20} \end{align}\] with \(H_{tc}\) given in Eq. 13 , \({}_{T_1}\xi_m\) and \({}_{T_2}\xi_m , \, m \in \{1,2,3\}\) the components of the complex values modified frequency vectors, given by \[{}_{T_1}\boldsymbol{\xi} = \begin{pmatrix} \displaystyle \frac{2}{h_1}\, \mathrm{sin}\left(\frac{\xi_1^{\mathrm{BCs}}}{2}\right)\,\mathrm{cos}\left(\frac{\xi_2^{\mathrm{BCs}}}{2}\right)\,\mathrm{cos}\left(\frac{\xi_3^{\mathrm{BCs}}}{2}\right) \\[1em] \displaystyle \frac{2}{h_2}\, \mathrm{cos}\left(\frac{\xi_1^{\mathrm{BCs}}}{2}\right)\,\mathrm{sin}\left(\frac{\xi_2^{\mathrm{BCs}}}{2}\right)\,\mathrm{cos}\left(\frac{\xi_3^{\mathrm{BCs}}}{2}\right) \\[1em] \displaystyle \frac{2}{h_3}\, \mathrm{cos}\left(\frac{\xi_1^{\mathrm{BCs}}}{2}\right)\,\mathrm{cos}\left(\frac{\xi_2^{\mathrm{BCs}}}{2}\right)\,\mathrm{sin}\left(\frac{\xi_3^{\mathrm{BCs}}}{2}\right) \end{pmatrix} - \mathrm{i} \begin{pmatrix} \displaystyle \frac{2}{h_1}\, \cos\left(\frac{\xi_1^{\mathrm{BCs}}}{2}\right)\,\sin\left(\frac{\xi_2^{\mathrm{BCs}}}{2}\right)\,\sin\left(\frac{\xi_3^{\mathrm{BCs}}}{2}\right) \\[1em] \displaystyle \frac{2}{h_2}\, \sin\left(\frac{\xi_1^{\mathrm{BCs}}}{2}\right)\,\cos\left(\frac{\xi_2^{\mathrm{BCs}}}{2}\right)\,\sin\left(\frac{\xi_3^{\mathrm{BCs}}}{2}\right) \\[1em] \displaystyle \frac{2}{h_3}\, \sin\left(\frac{\xi_1^{\mathrm{BCs}}}{2}\right)\,\sin\left(\frac{\xi_2^{\mathrm{BCs}}}{2}\right)\,\cos\left(\frac{\xi_3^{\mathrm{BCs}}}{2}\right) \end{pmatrix} \label{eq:freq95TETRA2}\tag{21}\] \[{}_{T_2}\boldsymbol{\xi} = \overline{{}_{T_1}\boldsymbol{\xi}}\] Here, the frequencies \(\xi_i^{\mathrm{BCs}}\), \(i \in \{1,2,3\}\) corresponds to the periodic case values. For a full periodic field, its gradient, divergence and Laplacian can be expressed. In this regard, for a node-defined displacement vector field \(\boldsymbol{v}\) and a voxel-defined second-order stress/strain tensor \(\oalign{\boldsymbol{S}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\), we can write \[\widehat{\oalign{\boldsymbol{{}_{T_r}\nabla v}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}} = \mathrm{i} \, \widehat{\boldsymbol{v}} \otimes {}_{T_r}\boldsymbol{\xi} \,H_{tc} \quad; \quad \widehat{\boldsymbol{{}_{T_r}\nabla}\cdot \oalign{\boldsymbol{S}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}} = \mathrm{i} \, \widehat{\oalign{\boldsymbol{S}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}}\,{}_{T_r}\boldsymbol{\xi} \,\overline{H_{tc}}\] \[\widehat{\boldsymbol{{}_{T_1T_2}\Delta v}} = \frac{1}{2}\left(\widehat{D_m^{T_1}D_m^{T_2}\,\boldsymbol{v}} + \widehat{D_m^{T_2}D_m^{T_1}\,\boldsymbol{v}}\right) = - \lVert {}_{T_1}\boldsymbol{\xi} \rVert^2 \, \widehat{\boldsymbol{v}} = - \lVert {}_{T_2}\boldsymbol{\xi} \rVert^2 \, \widehat{\boldsymbol{v}} \label{eq:laplacian95tetra2}\tag{22}\]

It is important to recall that for the second-order derivatives (e.g., Laplacian operator in Eq. 22 ), the issue associated with the first derivative do not hold (i.e., a field and its Laplacian possess the same symmetries). For fields with symmetry extensions, the second-order derivatives can then be evaluated using appropriate connection between DTTs and DFT on extended field (see Tab. 1 in Sec. 2.4). In this case, the frequencies \(\xi_i^{\mathrm{BCs}}\), \(i \in \{1,2,3\}\) used in Eq. 21 corresponds to the non-periodic values in Tab. 1.

2.6 Discrete Green operators↩︎

The previous sections elaborated on the accelerated fixed-point approach, as well as the evaluation of discrete derivatives (in real and Fourier space). The current section is dedicated to the expression of the discrete Green operator \({}_d\widehat{\mathrm{DGO}}\), i.e., in general words, the deduction of the fluctuation displacement field, \(\boldsymbol{u^f}\), knowing the divergence of polarization field, \(\boldsymbol{p}\) (see Eq. 6 ). In this work, we distinguish the discrete Green operator as function of the pre-conditioner \(\mathbb{R}_0\), the derivation scheme and the loading conditions. Classically, in the periodic case, the pre-conditioner \(\mathbb{R}_0\) (generally called “reference material”) is often chosen as an isotropic elastic stiffness tensor. In the non-periodic case, its usage is limited to specific BCs, and two alternatives [19], [20] are discussed.

2.6.1 \(\mathbb{C}_0\)-discrete Green operator↩︎

Within the finite transformation framework, we first consider the usual choice of pre-conditioner such as \[\mathbb{R}_0:{}_d\oalign{\boldsymbol{\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \mathbb{C}_0:{}_d\oalign{\boldsymbol{\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = 2\mu_0 \, {}_d\oalign{\boldsymbol{\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} + \lambda_0 \, \mathrm{tr}\left({}_d\oalign{\boldsymbol{\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\right) \boldsymbol{1} \label{eq:TF95C095behavior}\tag{23}\] with \(\mathbb{C}_0\) the stiffness-like matrix, \(\lambda_0\) and \(\mu_0\) respectively the first and second Lamé parameters. Choosing these parameters appropriately is crucial for ensuring convergence. Within the small-strain framework, Eq. 23 is updated by replacing the fluctuation displacement gradient \({}_d\oalign{\boldsymbol{\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\) with its symmetric part \({}_d\oalign{\boldsymbol{\nabla^s u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\). In that case, the subsequent analysis remain valid. Focusing on the finite transformation framework, the discrete Green operator, \({}_d\widehat{\mathrm{DGO}}\), is then deduced by detailing the following expression in Fourier space \[\begin{align} \widehat{\boldsymbol{p}} &= \widehat{\boldsymbol{{}_d\nabla}\cdot\left(\mathbb{C}_0: \oalign{\boldsymbol{{}_d\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\right)} \nonumber\\ & = 2\mu_0 \, \widehat{\boldsymbol{{}_d\nabla}\cdot\oalign{\boldsymbol{{}_d\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}} + \lambda_0 \, \widehat{\boldsymbol{{}_d\nabla}\cdot\oalign{\boldsymbol{{}_d\nabla^T u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}} \label{eq:FourierSGO} \end{align}\tag{24}\]

In the context of general (periodic and non-periodic) BCs, the calculation of \(\widehat{\boldsymbol{p}}\) can be done using DTs (depending on the corresponding symmetries and/or periodicity). Based on Eq. 24 , the periodicity and or symmetry extensions of \(\boldsymbol{p}\) are required to be the same as the ones of \(\boldsymbol{{}_d\nabla}\cdot\oalign{\boldsymbol{{}_d\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\) and \(\boldsymbol{{}_d\nabla}\cdot\oalign{\boldsymbol{{}_d\nabla^T u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\). For the sake of this investigation, we consider that a given component \(u^f_i,\,i \in \left\{1,2,3\right\}\) has the symmetries and/or periodicity \(A_iB_iC_i\), where \(A_i\), \(B_i\) and \(C_i\) are respectively the extensions in the three directions \(\boldsymbol{e}_1\), \(\boldsymbol{e}_2\) and \(\boldsymbol{e}_3\) (i.e., w.r.t. the faces \(S_{1\alpha}\), \(S_{2\alpha}\) and \(S_{3\alpha}\), \(\alpha \in \{0,1\}\)). Based on the previous analysis of derivative properties in Eq. 11 and examining the components \(\left[\boldsymbol{\nabla}\cdot\oalign{\boldsymbol{\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\right]_{j}=u^f_{j,ii}\) and \(\left[\boldsymbol{\nabla}\cdot\oalign{\boldsymbol{\nabla^T u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\right]_{j}=u^f_{i,ji}\), we can write the extensions of each term involved in these operations as follows \[\boldsymbol{\nabla}\cdot\oalign{\boldsymbol{\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \begin{pmatrix} \displaystyle \underbrace{u^f_{1,11}}_{A_1B_1C_1} + \underbrace{u^f_{1,22}}_{A_1B_1C_1} + \underbrace{u^f_{1,33}}_{A_1B_1C_1}\\[2em] \displaystyle \underbrace{u^f_{2,11}}_{A_2B_2C_2} + \underbrace{u^f_{2,22}}_{A_2B_2C_2} + \underbrace{u^f_{2,33}}_{A_2B_2C_2} \\[2em] \displaystyle \underbrace{u^f_{3,11}}_{A_3B_3C_3} + \underbrace{u^f_{3,22}}_{A_3B_3C_3} + \underbrace{u^f_{3,33}}_{A_3B_3C_3} \end{pmatrix} \quad; \quad \boldsymbol{\nabla}\cdot\oalign{\boldsymbol{\nabla^T u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \begin{pmatrix} \displaystyle \underbrace{u^f_{1,11}}_{A_1B_1C_1} + \underbrace{u^f_{2,12}}_{\overline{A_2B_2}C_2} + \underbrace{u^f_{3,13}}_{\overline{A_3}B_3\overline{C_3}}\\[2em] \displaystyle \underbrace{u^f_{1,21}}_{\overline{A_1B_1}C_1} + \underbrace{u^f_{2,22}}_{A_2B_2C_2} + \underbrace{u^f_{3,23}}_{A_3\overline{B_3C_3}} \\[2em] \displaystyle \underbrace{u^f_{1,31}}_{\overline{A_1}B_1\overline{C_1}} + \underbrace{u^f_{2,32}}_{A_2\overline{B_2C_2}} + \underbrace{u^f_{3,33}}_{A_3B_3C_3} \end{pmatrix} \label{eq:lapl95laplT95SGO95BCs}\tag{25}\] In Eq. 25 , the underbrace expression corresponds to the symmetries and/or periodicity of a given field. Based on the BC analysis in Eq. 25 , in the general case, no obvious symmetries and/or periodicity can be identified for the field \(\boldsymbol{\nabla}\cdot\oalign{\boldsymbol{\nabla^T u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\) and consequently on the field \(\boldsymbol{p}\). Nevertheless, choosing appropriately the BCs, it is still possible to constrain both terms \(\boldsymbol{\nabla}\cdot\oalign{\boldsymbol{\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\) and \(\boldsymbol{\nabla}\cdot\oalign{\boldsymbol{\nabla^T u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\) to have the same symmetries and/or periodicity. The appropriate choice to fulfill this condition is such as \[\left\{ \begin{array}{c} A_1 = \overline{A_2} = \overline{A_3} \\[0.2em] B_1 = \overline{B_2} = B_3 \\[0.2em] C_1 = C_2 = \overline{C_3} \end{array} \right. \label{eq:SGO95BCs}\tag{26}\] so that \(\boldsymbol{p}\) and \(\boldsymbol{u^f}\) exhibit identical symmetry and/or periodicity extensions. In practice, Eq. 26 imposes to the fluctuation displacement field \(\boldsymbol{u^f}\) to have the following periodicity and/or symmetry extensions:

1. per face, the full anti-symmetry of the normal components and full symmetry of the tangential ones,

2. or, conversely, per face, the full anti-symmetry of the tangential components and full the symmetry of the normal ones,

3. per opposite faces, the periodicity of the three components

These categories of conditions can be referred to the so-called “normal-mixed” BCs. In total the conditions in Eq. 26 lead to a maximum of \(5^3\) loading possibilities (5 possibilities per couple of opposite faces and per component). Loading conditions used in the work of [17] in the context of homogenization are particular cases of the “normal-mixed” BCs as introduced here. For “normal-mixed” BCs, depending on the finite difference scheme, the relation in Eq. 24 can be detailed.

• HEX1 scheme:

Considering the HEX1 scheme and combining the relations in Eq. 15 with Eq. 24 , we have \[-\widehat{\boldsymbol{p}} = 2\mu_0 \lVert {}_{H}\boldsymbol{\xi} \rVert^2 \widehat{\boldsymbol{u^f}} + \lambda_0 \left({}_{H}\boldsymbol{\xi} \otimes {}_{H}\boldsymbol{\xi}\right) \widehat{\boldsymbol{u^f}} = 2\mu_0 \lVert {}_{H}\boldsymbol{\xi} \rVert^2 \widehat{\boldsymbol{u^f}} + \lambda_0 \left({}_{H}\boldsymbol{\xi} \cdot \widehat{\boldsymbol{u^f}}\right) {}_{H}\boldsymbol{\xi}\] Multiplying this relation by \({}_{H}\boldsymbol{\xi}\) helps to get the term \(\widehat{\boldsymbol{u^f}}\cdot{}_{H}\boldsymbol{\xi}\) as a function of \(\widehat{\boldsymbol{p}}\cdot{}_{H}\boldsymbol{\xi}\). It follows that \[\widehat{\boldsymbol{u^{f}}} = - \displaystyle \frac{1}{2\mu_0 \lVert {}_{H}\boldsymbol{\xi} \rVert^2} \left(\widehat{\boldsymbol{p}} - \frac{ \lambda_0}{2\mu_0 + \lambda_0} \frac{\widehat{\boldsymbol{p}} \cdot {}_{H}\boldsymbol{\xi}}{ \lVert {}_{H}\boldsymbol{\xi} \rVert^2} {}_{H}\boldsymbol{\xi} \right) = {}_H\widehat{\mathrm{DGO}}(\widehat{\boldsymbol{p}}) \label{eq:FourierSGO95HEX1}\tag{27}\] This expression corresponds to the full periodic case [10], [13], with the notable distinction that the modified frequency vector \({}_{H}\boldsymbol{\xi}\) must account for the BCs as described in Eq. 13 .

• TETRA2 scheme:

With the TETRA2 scheme, due to the second-order derivatives in Eq. 24 , the DFT of the divergence of the polarization stress \(\widehat{\boldsymbol{p}}\) is written as \[\widehat{\boldsymbol{p}} = \displaystyle \frac{1}{2} \left(\widehat{{}_{T_2T_1}\boldsymbol{p}} + \widehat{{}_{T_1T_2}\boldsymbol{p}}\right)\] with \[\left\{ \begin{array}{l} \widehat{{}_{T_2T_1}\boldsymbol{p}} = 2\mu_0 \, \widehat{\boldsymbol{{}_{T_2}\nabla}\cdot\oalign{\boldsymbol{{}_{T_1}\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}} + \lambda_0 \, \widehat{\boldsymbol{{}_{T_2}\nabla}\cdot\oalign{\boldsymbol{{}_{T_1}\nabla^T u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}} \\[0.5em] \widehat{{}_{T_1T_2}\boldsymbol{p}} = 2\mu_0 \, \widehat{\boldsymbol{{}_{T_1}\nabla}\cdot\oalign{\boldsymbol{{}_{T_2}\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}} + \lambda_0 \, \widehat{\boldsymbol{{}_{T_1}\nabla}\cdot\oalign{\boldsymbol{{}_{T_2}\nabla^T u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}} \end{array} \right.\] Consequently, \[-\widehat{\boldsymbol{p}} = \displaystyle 2 \mu_0 \, \lVert {}_{T_2}\boldsymbol{\xi} \rVert^2 \, \widehat{\boldsymbol{u^f}} + \frac{\lambda_0}{2} \, \left[ \left(\widehat{\boldsymbol{u^f}}\cdot \overline{{}_{T_2}\boldsymbol{\xi}} \right) {}_{T_2}\boldsymbol{\xi} + \left(\widehat{\boldsymbol{u^f}}\cdot {}_{T_2}\boldsymbol{\xi} \right) \overline{{}_{T_2}\boldsymbol{\xi}} \right] \label{eq:p95u95tetra2}\tag{28}\] Multiplying Eq. 28 by \({}_{T_2}\boldsymbol{\xi}\) and \(\overline{{}_{T_2}\boldsymbol{\xi}}\) leads to the system of equations \[\left\{ \begin{array}{l} -\widehat{\boldsymbol{p}} \cdot {}_{T_2}\boldsymbol{\xi} = \displaystyle \frac{2\mu_0+\lambda_0}{2} \, \lVert {}_{T_2}\boldsymbol{\xi} \rVert^2 \, \left(\widehat{\boldsymbol{u^f}}\cdot {}_{T_2}\boldsymbol{\xi} \right) + \frac{\lambda_0}{2} \, \left({}_{T_2}\boldsymbol{\xi} \cdot {}_{T_2}\boldsymbol{\xi}\right) \, \left(\widehat{\boldsymbol{u^f}}\cdot \overline{{}_{T_2}\boldsymbol{\xi}} \right)\\[1em] -\widehat{\boldsymbol{p}} \cdot \overline{{}_{T_2}\boldsymbol{\xi}} = \displaystyle \frac{\lambda_0}{2} \, \left(\overline{{}_{T_2}\boldsymbol{\xi}} \cdot \overline{{}_{T_2}\boldsymbol{\xi}}\right) \, \left(\widehat{\boldsymbol{u^f}}\cdot {}_{T_2}\boldsymbol{\xi} \right) + \frac{2\mu_0+\lambda_0}{2} \, \lVert {}_{T_2}\boldsymbol{\xi} \rVert^2 \, \left(\widehat{\boldsymbol{u^f}}\cdot \overline{{}_{T_2}\boldsymbol{\xi}} \right) \end{array} \right.\] It is possible to invert this system by considering the terms \(\left(\widehat{\boldsymbol{u^f}}\cdot {}_{T_2}\boldsymbol{\xi} \right)\) and \(\left(\widehat{\boldsymbol{u^f}}\cdot \overline{{}_{T_2}\boldsymbol{\xi}} \right)\) as the unknowns. Introducing these quantities back in Eq. 28 , the Fourier discrete Green operator for TETRA2 is given by \[\widehat{\boldsymbol{u^{f}}} = \displaystyle - \frac{1}{2\mu_0\,\lVert {}_{T_2}\boldsymbol{\xi} \rVert^2} \left[\widehat{\boldsymbol{p}} + \displaystyle \frac{ \lambda_0}{2}\left(p_{\mathrm{int}1}\,\overline{{}_{T_2}\boldsymbol{\xi}} + p_{\mathrm{int}2}\,{}_{T_2}\boldsymbol{\xi} \right)\right] = {}_{T_1T_2}\widehat{\mathrm{DGO}}(\widehat{\boldsymbol{p}}) \label{eq:FourierSGO95TETRA2951}\tag{29}\] with \(p_{\mathrm{int}1}\) and \(p_{\mathrm{int}2}\) intermediate fields defined as \[\left\{ \begin{array}{l} p_{\mathrm{int}1} = \displaystyle \frac{1}{D_s} \left[\displaystyle -\frac{2\mu_0+\lambda_0}{2}\,\lVert {}_{T_2}\boldsymbol{\xi} \rVert^2\,\left( \widehat{\boldsymbol{p}} \cdot {}_{T_2}\boldsymbol{\xi}\right) + \displaystyle \frac{\lambda_0}{2}\, \left({}_{T_2}\boldsymbol{\xi} \cdot {}_{T_2}\boldsymbol{\xi}\right) \left( \widehat{\boldsymbol{p}} \cdot \overline{{}_{T_2}\boldsymbol{\xi}}\right)\right] \\[1em] p_{\mathrm{int}2} = \displaystyle \frac{1}{D_s} \left[ \displaystyle \frac{\lambda_0}{2}\,\left(\overline{{}_{T_2}\boldsymbol{\xi}} \cdot \overline{{}_{T_2}\boldsymbol{\xi}}\right)\,\left( \widehat{\boldsymbol{p}} \cdot {}_{T_2}\boldsymbol{\xi}\right) - \displaystyle\frac{2\mu_0+\lambda_0}{2} \, \lVert {}_{T_2}\boldsymbol{\xi} \rVert^2 \left( \widehat{\boldsymbol{p}} \cdot \overline{{}_{T_2}\boldsymbol{\xi}}\right)\right] \\[1em] \qquad \text{where } D_s = \displaystyle \left(\frac{2\mu_0+\lambda_0}{2}\,\lVert {}_{T_2}\boldsymbol{\xi} \rVert^2\right)^2 - \left(\frac{\lambda_0}{2}\,\lvert {}_{T_2}\boldsymbol{\xi} \cdot {}_{T_2}\boldsymbol{\xi} \rvert\right)^2 \end{array} \right. \label{eq:FourierSGO95TETRA2952}\tag{30}\]

5 shows the small strain version of the discrete Green operator in Eqs. 27 and 29 . It is important to recall that the discrete Green operators in Eqs. 27 and 29 are valid for “normal-mixed” BCs. However, it should be noted that, e.g., a bending loading BCs or a full Dirichlet BCs, are not included in this list of \(5^3\) eligible possibilities. To overcome this limitation, the reference material behavior must be defined alternatively. Two options are examined hereafter.

2.6.2 \(2\mu_0\)-discrete Green operator↩︎

As discussed in the previous section, the restriction to specific BCs (see Eq. 26 ) stems from the symmetries and/or periodicity of certain derivative terms appearing in the computation of \(\widehat{\boldsymbol{\nabla}\cdot\oalign{\boldsymbol{\nabla^T u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}}\) in Eq. 24 . To overcome this limitation, a simple choice is to set the first Lamé parameter to zero, i.e., \(\lambda_0=0\) [19], [20]. This means that the pre-conditioner is described by \[\mathbb{R}_0:\oalign{\boldsymbol{{}_d\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = 2\mu_0\,\oalign{\boldsymbol{{}_d\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\] This choice eliminates the need for any constraints on the periodicity or symmetry extensions of the fluctuation displacement field \(\boldsymbol{u^f}\). The result of this process is the possible combination of any Neumann, Dirichlet and periodic BCs for all faces and components, which yields a total of \(5^9\) possibilities.

Exploiting the previous analysis and fixing \(\lambda_0=0\), we deduce the Fourier discrete Green operator as

• HEX1 scheme case \[\widehat{\boldsymbol{u^{f}}} = - \displaystyle \frac{1}{2\mu_0 \lVert \boldsymbol{\xi}_{H} \rVert^2} \widehat{\boldsymbol{p}} = {}_{H}\widehat{\mathrm{DGO}}(\widehat{\boldsymbol{p}}) \label{eq:NSGO195HEX1}\tag{31}\]

• TETRA2 scheme case \[\widehat{\boldsymbol{u^{f}}} = - \displaystyle \frac{1}{2\mu_0 \lVert {}_{T_2}\boldsymbol{\xi} \rVert^2} \widehat{\boldsymbol{p}} = {}_{T_1T_2}\widehat{\mathrm{DGO}}(\widehat{\boldsymbol{p}}) \label{eq:NSGO195TETRA2}\tag{32}\]

The constructed discrete Green operators are analogous for both finite difference schemes. This formulation offers the advantage of simplicity, requiring only the choice of the second Lamé parameter, \(\mu_0\), and is suitable for all types of BCs.

2.6.3 \(\mathbb{B}_0\)-discrete Green operator↩︎

In contrast with the previous approach of setting the first Lamé parameter of the reference material to zero, another potential solution is investigated in this section. Following the work of [20], an alternative pre-conditioner is \[\mathbb{R}_0: \oalign{\boldsymbol{\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = 2\mu_0 \oalign{\boldsymbol{\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} + \lambda_0 \oalign{\boldsymbol{\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} \odot \oalign{\boldsymbol{I}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} =\mathbb{B}_0: \oalign{\boldsymbol{\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\] Similarly to the previous cases, the discrete Green operator \({}_d\widehat{\mathrm{DGO}}\), for a given finite difference scheme, is determined based on the following relationship \[\widehat{\boldsymbol{p}} = \widehat{\boldsymbol{{}_d\nabla}\cdot\left(\mathbb{B}_0: \oalign{\boldsymbol{{}_d\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\right)} = 2\mu_0 \, \widehat{\boldsymbol{{}_d\nabla}\cdot\oalign{\boldsymbol{{}_d\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}} + \lambda_0 \, \widehat{\boldsymbol{{}_d\nabla}\cdot\left(\oalign{\boldsymbol{{}_d\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} \odot \oalign{\boldsymbol{I}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\right)} \label{eq:FourierNSGO2}\tag{33}\] In order to detail this discrete Green operator, the index formulation of Eq. 33 is used \[\begin{align} \widehat{p_i} &= \displaystyle 2 \mu_0 \, \sum_{j=1}^{3} \widehat{D_j^d D_j^d u^f_i} + \lambda_0 \, \sum_{j=1}^{3} \widehat{D_j^d\left(D_j^d u^f_i \, \delta_{ij}\right)} \qquad \text{ (no summation on the index } i \text{)}\\ & = \displaystyle 2 \mu_0 \, \sum_{j=1}^{3} \widehat{D_j^d D_j^d u^f_i} + \lambda_0 \, \widehat{D_i^d D_i^d u^f_i} \qquad \text{ (no summation on the index } i \text{)} \end{align}\] This expression indicates that the term \(\boldsymbol{{}_d\nabla}\cdot\left(\oalign{\boldsymbol{{}_d\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} \odot \oalign{\boldsymbol{I}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\right)\) yields equivalent periodicity and/or symmetry extensions to those of \(\boldsymbol{{}_d\nabla}\cdot\oalign{\boldsymbol{{}_d\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\), which are also identical to the ones of \(\boldsymbol{u^f}\). Consequently, the constructed discrete Green operator can be used for any choice of BCs.

• Accordingly to the HEX1 scheme and building on Eq. 12 , we have

\[\widehat{D_j^H D_j^H u^f_i} = - \lVert {}_{H}\boldsymbol{\xi} \rVert^2 \, \widehat{u^f_i}\] \[\widehat{D_i^H D_i^H u^f_i} = - \left({}_H \xi_i\right)^2 \, \widehat{u^f_i} \qquad \text{ (no summation on the index } i \text{)}\] It follows that \[\widehat{p_i} = - 2 \mu_0 \, \lVert {}_{H}\boldsymbol{\xi} \rVert^2 \, \widehat{u^f_i} - \lambda_0 \, \left({}_H \xi_i\right)^2 \, \widehat{u^f_i} \qquad \text{ (no summation on the index } i \text{)}\] Finally, \[\widehat{u^f_i} = - \displaystyle \frac{1}{2\mu_0 \lVert {}_{H}\boldsymbol{\xi} \rVert^2 + \lambda_0 \, \left({{}_{H}\xi}_i\right)^2}\, \widehat{p_i} = {}_{H}\widehat{\mathrm{DGO}}(\widehat{\boldsymbol{p}}) \qquad \text{ (no summation on the index } i \text{)} \label{eq:NSGO295HEX1}\tag{34}\]

• Using the same previous calculation steps for the TETRA2 scheme, we obtain the discrete Green operator as

\[\widehat{u^f_i} = - \displaystyle \frac{1}{2\mu_0 \lVert {}_{T_2}\boldsymbol{\xi} \rVert^2 + \lambda_0 \, \lvert {{}_{T_2}\xi}_i\rvert^2} \,\widehat{p_i} = {}_{T_1T_2}\widehat{\mathrm{DGO}}(\widehat{\boldsymbol{p}}) \qquad \text{ (no summation on the index } i \text{)} \label{eq:NSGO295TETRA2}\tag{35}\]

Unless otherwise specified, the discrete Green operator based on the pre-conditioner, \(\mathbb{B}_0: \oalign{\boldsymbol{\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\) (reference material) is the one used in the simulations (small and finite transformations). It is worth noting that for all the three versions of the discrete Green operators (see Eqs. 27 , 29 30 , 31 , 32 , 34 and 35 ), when the corresponding denominators are found to be null, we enforce \(\widehat{\boldsymbol{u^{f}}} = 0\).

2.7 Implementation specificities↩︎

The purpose of this section is to emphasize the most important evolutions (i.e., w.r.t. a classical periodic implementation) required to introduce non-periodic BCs together with various types of finite difference schemes in a distributed memory parallel solver.

2.7.1 Field implementation↩︎

• Fields in real space

In most classical implementations of periodic FFT-based solvers, fields in real space are defined as arrays of size (\(n_1^v\),\(n_2^v\),\(n_3^v\),\(n_{comp}\)), with \(n_{comp}\) the number of components (3 for a vector, 6 for a symmetric second-order tensor and 9 for a non-symmetric tensor). The displacement and strain/stress fields are defined on grids \(\llbracket 0,n_1^v-1\rrbracket \times \llbracket 0,n_2^v-1\rrbracket \times \llbracket 0,n_3^v-1\rrbracket\).

In our implementation for generic BCs, the strain field is still defined on a grid \(\llbracket 0,n_1^v-1\rrbracket \times \llbracket 0,n_2^v-1\rrbracket \times \llbracket 0,n_3^v-1\rrbracket\) while the displacement is now defined on a larger grid \(\llbracket 0,n_1^v\rrbracket \times \llbracket 0,n_2^v\rrbracket \times \llbracket 0,n_3^v\rrbracket\). In addition, for an extension to various types of finite difference schemes [35] requiring more than one point per voxel, fields are defined as arrays of size (\(n_1^v\),\(n_2^v\),\(n_3^v\),\(n_{comp}\),\(n_{ppv}\)), with \(n_{ppv}\) the number of points per voxel (e.g., using TETRA2 scheme, \(n_{ppv}=2\) for a strain/stress field and 1 for a displacement field). As a consequence, new field objects are defined to distinguish the two kinds of fields and introduce multi-points per voxel.

• Fields in Fourier space

In classical (periodic) implementations, all fields in Fourier space are complex numbers with a size (\(n_1^v/2+1\),\(n_2^v\),\(n_3^v\),\(n_{comp}\)). Actually, to reduce memory footprint and building on the fact that the input fields are real values, the size of the DFT in the first direction is classically divided by (almost) two without any loss of information (only non-redundant values are stored).

To account for possibly non-periodic BCs, a field component in Fourier space may be obtained either through a DFT, when periodicity is prescribed in at least one direction, or through DTTs when non-periodic BCs are considered. Consequently, depending on the selected BCs, a component of a field, in Fourier space, may be either complex-valued (if a DFT is involved in at least one direction) or purely real-valued (if DTTs are employed in all three directions). In addition, if a periodicity condition appears firstly in a given direction \(l\in \{1,2,3\}\), the array dimension divided by \(2\) to save memory corresponds to the precise direction \(l\). As a consequence, to propose a versatile code accounting for general BCs (i.e., arbitrary combinations of periodic and/or non-periodic BCs on each face and each component), the novel implementation in this work includes Fourier field object that are able to handle components of different types (complex or real) and sizes.

2.7.2 Domain decomposition and discrete transformations↩︎

The distributed memory implementation of the FFTW library [29] is based on a 1D slab decomposition of the unit-cell. Such an implementation exhibits an important limit regarding the maximum number of processors. Actually, a \((n^v)^3\) unit-cell can be distributed on a maximum of \(n^v\) processors. To overcome this limitation, the 2DECOMP&FFT library [30], [31] provides the possibility of using 2D pencils to decompose the unit-cell. Fig. 6 shows an example of decomposition with \(12\) processors, each processor managing a part of the unit-cell. This MPI-based library is a foundation of the classical periodic solver [28], and the present work necessitated an evolution of 2DECOMP&FFT to account for DTTs in a massively parallel codes. In AMITEX\(^\star\), for field components in real space, the pencil direction is \(\boldsymbol{e}_1\). The 3D DTs consist of a succession of 1D transforms in each direction. The first transform is in the direction \(\boldsymbol{e}_1\). The second transform, in direction \(\boldsymbol{e}_2\), can not be made directly because \(e_1\)-pencils do not have the full length of data in direction \(\boldsymbol{e}_2\). Hence, a transposition of the array is done (via a MPI_ALLTOALL communication), before applying the DTs, to get it distributed in \(e_2\)-pencils (see Fig. 6). After the application of the DT in the third direction, requiring a second transposition step, the array is distributed in \(e_3\)-pencils. Consequently, field components in Fourier space are distributed in \(e_3\)-pencils.

Recalling that the displacement field is defined on the grid \(\llbracket 0,n_1^v\rrbracket \times \llbracket 0,n_2^v\rrbracket \times \llbracket 0,n_3^v\rrbracket\), it is important to notice that, within the novel AMITEX\(^\star\) solver, the application of the DTs in each direction must be done taking into account “dismissed” points (see index \(j\) Tab. 1 in Sec. 2.4). For example, consider the sine transform DST1 of a 1D signal, the value of the signal at the two sides are set to zero and the DST1 is applied to a reduced signal (skipping the two extreme points). Depending on the DT, the first and last point could be dismissed or not. The “dismissed” points are identified during AMITEX\(^\star\) initialization and provided to 2DECOMP&FFT when initializing a 3D transform in the new library.

Actually, when a 2DECOMP&FFT 3D transform is initialized, an array of integers is provided to describe the DTs in the three directions. The first three components of the array are mandatory and describe the forward DTs in the three directions \(\boldsymbol{e}_1\), \(\boldsymbol{e}_2\) and \(\boldsymbol{e}_3\). Additionally, optional arguments can be provided to describe the number of elements “dismissed” in each direction, and their respective location (left and/or right side). The available 1D transforms are DFT, alongside with 8 DTTs (for the present displacement-based algorithm, only \(4\) DTTs are used), leading to \(9^3 = 729\) combinations for the 3D transform. The DTTs available in the 2DECOMP&FFT library correspond exactly to the transforms FFTW_REDFT (i.e., cosine transforms) and FFTW_RODFT (i.e., sine transforms) available in FFTW3 [29].

The internal memory requirement for the parallel 3D transform in 2DECOMP&FFT can reach five 3D temporary arrays. A reduction of this memory usage is possible, but it would increase the complexity of the 2DECOMP&FFT code. To limit the impact of this increased memory usage, the library now provides a memory pool. When the caller (i.e., AMITEX\(^\star\) in the present case) is not performing 3D transforms, he can temporarily read and write the memory blocks allocated by 2DECOMP&FFT. Note that, the library being open source, all these developments are available to anyone, especially for applications in fluid/solid mechanics and multi-physics.

The implementation of the DTTs done in the parallel framework (MPI), have been validated by the following two-step methodology. First, considering a random 3D field (scalar \(G\), vector \(\boldsymbol{G}\) or tensor \(\oalign{\boldsymbol{G}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\)) with various symmetries, the corresponding DTT computation of the field is undertaken to obtain a field in spectral space. Next, an inverse DTT calculation, denoted iDTT, is performed on the obtained spectral field. This sequence of operations (e.g., in the vector case iDTT(DTT)(\(\boldsymbol{G}\))) is expected to result in the recovery of the original random field (e.g., in the vector case \(\boldsymbol{G}\)). To ensure completeness, this validation is performed for any combination (i.e., symmetries and/or periodicity) of DTs (one per direction) varying also the grid sizes and the domain decomposition (number of processors).

2.7.3 Discrete derivation↩︎

• Derivation in real space

The present displacement-based algorithm (see Alg. 4) differs from the classical strain-based algorithm used in the full periodic case [2], [28]. Indeed, the displacement-based algorithm has the advantage to reduce the memory footprint (by a factor of \(2\) or \(3\) whether small or finite transformation hypothesis are considered) of the array used to store the \(4\) couples of solutions and residual fields used for the Anderson convergence acceleration. A further distinction between the displacement-based algorithm and the strain-based version implemented in [28], is that the former involves the evaluation of the divergence and the gradient in real space (see Alg. 4), whereas these derivation operations are not computed in the latter.

In the present distributed memory implementation, using pencil decomposition (see previously described Fig. 6), computation of discrete derivation at points located at the frontier of the pencil domains can not be directly done. Indeed, the data of points located in the neighboring pencil are required for this derivation operation. Hence, such discrete derivation necessitates the use of halo-voxels: each pencil is then enlarged by one voxel (one is enough for HEX1 and TETRA2 schemes) and the corresponding data are transferred from one processor to another through a 2DECOMP&FFT procedure relying on MPI non-blocking send/receive functions. In addition, the halo-voxels that are defined outside of the unit-cell must be evaluated according to the appropriate symmetry extension (associated with BCs). This last point is performed at AMITEX\(^\star\) side.

• Derivation in Fourier space

The application of discrete derivation in Fourier space using DTs (see Sec. 2.5) requires the knowledge of modified frequency vectors defined in Eqs. 13 , 21 and Tab. 1 (that exhibits the links between DTTs and DFT). These modified frequency vectors are Fourier fields which are evaluated and stored “on the fly” as soon as required (for example when applying discrete divergence, gradient, or Green operators). Please note that, for the present algorithm, a modified frequency vector is stored per combination of BCs applied to the displacement field components. For full Dirichlet (or Neumann, or periodic) BC, the \(3\) components of the displacement have the same BCs and a single vector is stored. In the worst case, the three components of the displacement have different BCs and three modified frequency vectors are stored.

• Cross validation

As a result of the extension to non-periodic BCs, the implementation becomes more complex, prone to bugs, especially when considering the high number of possible BC combinations (i.e., the total number of BC combinations is \((5^3)^{n_{comp}}\), with \(n_{comp}=3\) for a vector displacement field, hence approximately \(2.10^6\) combinations). The validation of discrete derivation in real and Fourier spaces is essential. For that purpose, a cross validation strategy is implemented: a random field is initialized, then derived in real space on the one hand, and derived using DTs for the back and forth in Fourier space on the other hand. The relative squared error between the two derived fields is subsequently computed and verified to be lower than a tolerance of approximately \(10^{-13}\). This strategy is repeated for all the possible combinations of BCs, and for, scalar, vector and tensor (symmetric or not) as inputs of the various discrete derivation operators (divergence, gradient, curl, Laplacian). We believe that this cross validation is a major cornerstone in the development of the code: no need to go further if such cross validations are not performed.

3 Numerical simulations↩︎

A versatile parallel FFT-based solver, relying on the use of DTs, was introduced in the previous section to circumvent the limitations of periodic BCs. In the present section, an extensive investigation is conducted to highlight the extended capabilities and the robustness of this new solver. The objective is to demonstrate the versatility of the approach in addressing small and finite transformation problems with different loading scenarios of interest in material science and out of reach for standard (periodic) FFT-based solvers. The finite difference schemes introduced in this study are examined, and their outcomes are discussed, with a focus on the recently proposed TETRA2 scheme [27], that demonstrates a better robustness than the HEX1 scheme [24]. Leveraging on the parallel implementation, different types of unit-cell geometries and material behavior laws, including isotropic elasticity, isotropic perfect plasticity and crystal plasticity are explored. To ensure comprehensiveness, the behavior laws for each category of examples are briefly outlined in the corresponding section. When possible, the simulations are validated against analytical solutions. Unless otherwise precised, a tolerance \(\epsilon = 10^{-5}\) is used for the equilibrium condition (see Eq. 7 ).

3.1 Elastic cantilever beam under bending configurations↩︎

With the aim of validating non-homogeneous loadings, the first category of loadings, consists of a bending loading (see 6 for a homogeneous pure tension displacement-controlled analytical validation within finite transformation framework). Small and finite transformation frameworks are examined here using an elastic cantilever beam under two types of bending configurations. The small strain case for both the HEX1 and TETRA2 finite difference schemes is validated through the use of an analytical solution. The objective of the finite transformation case is to ascertain the most suitable candidate between the two presented finite difference schemes for such simulations. For simplicity, in the pure isotropic elastic cases, we consider standard steel elastic properties with the parameters: \(\lambda = 120 \, \si{GPa}\) and \(\mu = 80 \, \si{GPa}\). As the bending configurations to be treated are not compatible with “normal-mixed” (see Eq. 26 ), we then use the \(\mathbb{B}_0\)-discrete Green operator with \(\lambda_0=\lambda\) and \(\mu_0=\mu\).

3.1.1 Cantilever beam under uniformly distributed applied stress↩︎

Assuming small strain hypothesis and isotropic elasticity, the material behavior is given by \[\oalign{\boldsymbol{\sigma}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \mathbb{C}:\oalign{\boldsymbol{\nabla^s u}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \mu \, \oalign{\boldsymbol{\nabla^s u}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} + \lambda \, \mathrm{tr}(\oalign{\boldsymbol{\nabla u}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}) \, \oalign{\boldsymbol{1}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\] The well-known bending setup in Fig. 7 (a) is used. The left side (face \(S_{10}\)) of the beam is fully clamped by enforcing Dirichlet BCs with a zero value for the applied displacement field, i.e., \(\boldsymbol{u^*}=0\). On the top face \(S_{21}\), a Neumann BC with an uniform applied stress vector in the direction \(\boldsymbol{e}_2\) is imposed (i.e., \(\boldsymbol{t^*} = -t_2^* \boldsymbol{e}_2\)).

Considering that the length \(L\) of the beam in the direction \(\boldsymbol{e}_1\), is considerably greater than the other dimensions (\(h\) and \(w\) in the directions \(\boldsymbol{e}_2\) and \(\boldsymbol{e}_3\) respectively, thus allowing for the face load to be approximated by a linear load), the beam theory provides the analytical displacement solution, in the direction \(\boldsymbol{e}_2\), of the neutral axis \[u_2(x) = - \displaystyle \frac{t_2^*w}{EI} \left(\frac{L^2 x^2}{4} - \frac{Lx^3}{6} + \frac{x^4}{24}\right) \label{eq:analytical95sol95bending95elastic95own95weight}\tag{36}\] where \(x\) denotes here the position along \(\boldsymbol{e}_1\), \(I = wh^3/12\) is the area moment of inertia and \(E\) is the Young modulus. Fig. 7 (b) presents the comparison between the analytical solution (see Eq. 36 )and the FFT-based solver results for both finite difference schemes HEX1 and TETRA2. Three increasing loading steps are used to show that, within small strain framework, HEX1 and TETRA2 schemes give similar results in a very good agreement with the analytical solution.

3.1.2 Cantilever beam under end load applied stress↩︎

To check the robustness of the method, another kind of bending loading is investigated. More specifically, in addition, to small strain simulations, we also explore the capability of the solver for finite transformation simulations. The material behavior law is extended to a finite transformation formalism with the second Piola–Kirchhoff stress tensor given by the relationship \[\oalign{\boldsymbol{\Pi}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \mathbb{C}:\oalign{\boldsymbol{E}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \mu \, \oalign{\boldsymbol{E}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} + \lambda \, \mathrm{tr}(\oalign{\boldsymbol{E}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}) \, \oalign{\boldsymbol{1}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\] with the Green-Lagrange strain tensor \(\oalign{\boldsymbol{E}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\) defined as \[\oalign{\boldsymbol{E}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \displaystyle \frac{1}{2} \left(\oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^{T}\oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} - \oalign{\boldsymbol{1}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\right)\]

For the purpose of the present study, the bending setup with an end load applied stress as sketched in Fig. 8 (a) is used. This configuration prompts to a large displacement response as the loading increases. Fig. 8 (b) presents the macroscopic responses. In order to validate the implemented Lagrangian formulation for the finite transformation case, a comparison with the small strain framework simulations is also plotted. Fig. 8 (b) (left part) shows that both small and finite transformations lead to the same behavior at low applied stress. As expected, the choice of a small strain approximation leads to linear response (identical for TETRA2 and HEX1 schemes). Conversely, in the finite transformation framework, non-linear elastic responses are obtained. Both HEX1 and TETRA2 schemes yield similar non-linear behavior, allowing to reach more than \(0.5\) normalized displacement. It is important to notice that the finite transformation simulation with HEX1 scheme fails to converge post a certain level of applied stress \(T_{2\,\mathrm{failed}}^*\). Meanwhile, the TETRA2 scheme allows to converge well beyond the stress level \(T_{2\,\mathrm{failed}}^*\) (up to a ratio \(u_2(L)/L > 0.8\)). The total displacement profiles are plotted for both schemes at \(T_{2\,\mathrm{failed}}^*\), c.f., Fig. 8 (b) (right part). It is visible that the HEX1 scheme shows spurious displacement oscillations, leading to the failed convergence while TETRA2 scheme yields smooth deformed configuration, enabling convergence at elevated strain levels. The TETRA2 scheme seems to be more robust than the HEX1 scheme in modeling such non-linearity.

3.2 Perfect plasticity porous medium under a normal-mixed loading↩︎

Following the examples in elasticity, a perfect isotropic plasticity behavior in a finite transformation framework is considered in this section in order to study the deformation of a porous medium (see Fig. 9). FFT-based simulations of porous media in full periodic conditions have been reported in the literature [36], [37] using [28], to reproduce cavity growth and coalescence. These kind of numerical simulations are quite challenging. We discuss here, in a non-periodic case, the effect of the finite difference scheme, of the discrete Green operator and of the spatial discretization. For the simulations in this section only, a tolerance of \(\epsilon = 10^{-4}\) is employed for the equilibrium condition.

For the purpose, the small strain version of the classical von Mises isotropic \(J_2\) plasticity formulation is extended to finite transformation based on the logarithmic formalism [38]. This requires the introduction of the Hencky strain tensor as \[\oalign{\boldsymbol{H}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \frac{1}{2} \log \left(\oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^{T}\oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\right)\] Similarly to small strain framework an additive decomposition into an elastic part \(\oalign{\boldsymbol{H}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^e\) and a plastic contribution \(\oalign{\boldsymbol{H}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^p\) is considered for the Hencky strain tensor \[\oalign{\boldsymbol{H}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \oalign{\boldsymbol{H}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^e + \oalign{\boldsymbol{H}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^p\] The Cauchy stress is written as \[\oalign{\boldsymbol{\sigma}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = 2\mu \, \mathrm{dev}\left(\oalign{\boldsymbol{H}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^e\right) + K \, \mathrm{tr}(\oalign{\boldsymbol{H}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^e) \, \oalign{\boldsymbol{1}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\] with \(K = \lambda + 2 \mu /3\). Perfect plasticity is accounted for through a yield surface defined as \[f(\oalign{\boldsymbol{\sigma}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}) = \sigma_{eq} - \sigma_0 = 0\] where \(\sigma_{eq}\) is the equivalent stress written as \[\sigma_{eq} = \displaystyle \sqrt{\frac{3}{2} \, \oalign{\boldsymbol{s}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}:\oalign{\boldsymbol{s}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}} \qquad ; \qquad \oalign{\boldsymbol{s}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \mathrm{dev}\left(\oalign{\boldsymbol{\sigma}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\right)\] The normality flow rule is given by \[\dot{\oalign{\boldsymbol{H}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^p} = \dot{p} \dfrac{\partial f}{\partial \oalign{\boldsymbol{\sigma}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}} \qquad ; \qquad \dot{p} \geq 0\] with \(\dot{p}\) the rate of the cumulative plastic strain. The behavior law is implemented with MFront [39]. We use Young modulus \(E=200 \, \mathrm{GPa}\), a Poisson ratio \(\nu=0.3\) (i.e., \(\lambda\approx 115 \, \mathrm{GPa}\), \(\mu \approx 76 \, \mathrm{GPa}\)) and the initial yield stress \(\sigma_0 = 200\,\mathrm{MPa}\). The void is represented by an elastic material with null Lamé coefficients. The pre-conditioner parameters \(\lambda_0\) and \(\mu_0\) have to be chosen. To the best knowledge of the authors, no obvious optimized values for \(\lambda_0\) and \(\mu_0\) are demonstrated in the literature for non-periodic BCs in plasticity and finite transformation. Classical Moulinec and Suquet [2] choice would result in \(\beta_0 = 1/2 \left(\max (\beta) + \min (\beta)\right)\) with \(\beta \in \{\lambda,\mu\}\). Nevertheless, for the present simulations, based on preliminary tests (not shown here), we choose the fixed value of \(\beta_0=\max \left(\beta\right)\) (i.e., \(\lambda_0 = \max (\lambda)\) and \(\mu_0 = \max (\mu)\)) for all three pre-conditioners.

In the present non-periodic case, a homogenized-like tension which satisfies the “normal-mixed” condition in Eq. 26 is used. This gives the possibility to employ all three types of discrete Green operators. The BCs consist of normal-mixed BCs of type 1 (NMBC1) with Dirichlet conditions on the normal displacement components and Neumann conditions of the tangential stress, i.e., \[\text{NMBC1} \,\, \left\{ \begin{array}{l} \bullet \text{ Dirichlet BCs in directions } q=i, \quad u^f_i = 0 \text{ and } \\[0.5em] \bullet \text{ Null Neumann BCs in directions } q \neq i, \quad P_{iq} = 0 \end{array} \right. \label{eq:NMBC1}\tag{37}\] The loading is prescribed by imposing a volume-average of the gradient of displacement \(\oalign{\boldsymbol{\nabla u^*}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \oalign{\boldsymbol{G}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\) with the component \(G_{11}\) controlled and the other components \(G_{rs}\), \(\{r,s\} \neq \{1,1\}\) adjusted at each iteration to obtain \(\left<P_{rs}\right> = 0 \text{ if } \{r,s\} \neq \{1,1\}\).

Considering both finite difference schemes (HEX1 and TETRA2) and also various spatial discretizations, Fig. 10 gives the evolution of the average stress component, \(\left<P_{11}\right>\), as a function of the average of the displacement gradient component, \(\left<F_{11}-1\right>\). Overall, three regimes can be distinguished. The first regime consists of an elastic evolution up to an apparent yield stress \(\sigma_y\). It is followed by a softening regime due to void growth with a significant irreversible deformation. A quasi-linear decrease of the stress field from \(\sigma_y\) to a macroscopic “coalescence-like” stress \(\sigma_{c}\) is observed. The third regime characterizes an instability due to a strong plastic strain localization around the void, with a sudden drop of the macroscopic stress post-\(\sigma_{c}\). Going forward with the HEX1 scheme, increasing the discretization (from \(32^3\) to \(512^3\) voxels) leads to a visible increase in the apparent yield stress \(\sigma_y\) (see the zoom part in Fig. 10). On the contrary with the TETRA2 scheme, the apparent yield stress \(\sigma_y\) varies less with the discretization. Continuing with Fig. 10, we can see that both TETRA2 and HEX1 schemes, predict identical level of “coalescence-like” stress \(\sigma_{c}\) regardless of the discretization. It is possible to conclude that, with a reasonable discretization (e.g., \(64^3\) voxels), acceptable approximations of \(\sigma_y\) and \(\sigma_{c}\) are obtained with the TETRA2 scheme. Furthermore, examining the deformation level corresponding to the beginning of the “coalescence-like” transition, the TETRA2 scheme shows a higher sensitivity than the HEX1 scheme w.r.t. the discretization. However, with an increasing discretization, both the TETRA2 and HEX1 schemes seem to converge towards a macroscopic response. The TETRA2 scheme has the advantage of allowing convergence up to larger deformations in the “coalescence-like” regime in comparison with the HEX1 scheme.

Fig. 11 presents for a given discretization, i.e., \(n_x\times n_y \times n_z = 512^3\) voxels, the comparison of the three presented pre-conditioners (leading to the choice of discrete Green operator) considering the HEX1 and TETRA2 finite difference schemes. Globally, the invariance of the macroscopic responses is demonstrated with the selection of the pre-conditioners. More precisely, the simulations with the pre-conditioner \(2 \mu_0\oalign{\boldsymbol{\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\) failed to converge first, then follows the ones with the pre-conditioner \(\mathbb{C}_0: \oalign{\boldsymbol{\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\) and finally the simulations with the pre-conditioner \(\mathbb{B}_0: \oalign{\boldsymbol{\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\). The latter pre-conditioner helps to reach higher level of deformations. Tab. 2 summarizes the consequence of the choice of a discrete Green operator based on the finite difference scheme. Focusing on a deformation range where all three discrete Green operators converge, we can say that, for a given pre-conditioner, in order to reach a desired deformation, the TETRA2 scheme converges faster than the HEX1 scheme (approximatively \(1.4\) factor in terms of the iteration numbers). Based on a given finite difference scheme, the pre-conditioners \(\mathbb{C}_0: \oalign{\boldsymbol{\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\), \(\mathbb{B}_0: \oalign{\boldsymbol{\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\) and \(2 \mu_0\oalign{\boldsymbol{\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\) converge respectively one faster than the other. We recall that the pre-conditioner \(\mathbb{C}_0: \oalign{\boldsymbol{\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\) can only be used for BCs that are compatible with “normal-mixed” BCs in Eq. 26 .

To sum up, these results demonstrate the robustness of the novel FFT-based solver in handling complex simulations (non-periodic, infinite contrast, perfect plasticity, finite transformation). The TETRA2 scheme seems to demonstrate greater robustness than HEX1, in predicting the apparent yield stress (at small strain). However, the TETRA scheme exhibits higher discretization sensitivity in the “coalescence-like” strain level than the HEX1 scheme. A deeper analysis of the performance of the TETRA2 scheme and an optimization of the parameters \(\lambda_0\) and \(\mu_0\) are of interests but are beyond the scope of the present paper. It would also be interesting in a future work to enhance this study by including multiple voids with different shapes and apply the composite voxels technique [40], [41] specially in TETRA2 case.

3.3 Single crystal L-beam under torsion-bending loading↩︎

As it is thoroughly documented in the literature, the classical FFT-based solver has proven to be a valuable modeling technique for periodic single crystal and polycrystal materials [42], [43]. It is worthwhile to investigate the behavior of such anisotropic materials within the novel non-periodic framework. We aim to discuss observations within continuum single crystal plasticity, building on experimental studies such as shear-compression [44], torsion-bending [45], non-proportional orthogonal bending [46]. As part of the present work, these complex experimental-like simulations have been successfully tested. For conciseness, only the L-beam torsion-bending case inspired from the work of [45], is reported in this paper. Note that, we do not attempt here a direct point-to-point comparison with the experimental macroscopic curves. Such a comparison would require the utilization of higher-order type models (see e.g., a strain gradient crystal plasticity [47]–[49] or a Cosserat continuum framework [50]). These types of gradient-based plasticity theories are out of the scope of this paper and will be the subject of future work. A local crystal viscoplasticity model, which is based on the evolution of dislocation densities [51], [52], is employed for the purpose of the present analysis.

The finite transformation crystal plasticity model used in [52], is briefly recalled here. The deformation gradient, \(\oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \oalign{\boldsymbol{1}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} + \oalign{\boldsymbol{\nabla u}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\), is decomposed in a product of an elastic and a plastic parts [53], respectively noted \(\oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^e\) and \(\oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^p\) \[\oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^{e}\,\oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^{p}\] so that the Green-Lagrange elastic strain measure \(\oalign{\boldsymbol{E}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^e\) is written as \[\oalign{\boldsymbol{E}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^e = \displaystyle \frac{1}{2} \left(\oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^{eT}\oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^e - \oalign{\boldsymbol{1}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\right)\] In order to account for plasticity on each slip system, the total velocity gradient \(\oalign{\boldsymbol{L}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}=\dot{\oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}} \oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^{-1}\) is additively decomposed as \[\oalign{\boldsymbol{L}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \oalign{\boldsymbol{L}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^e + \oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^e \oalign{\boldsymbol{L}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^p {\oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^e}^{-1}\] with \[\oalign{\boldsymbol{L}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^{e} = \dot{\oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^{e}} {\oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^{e}}^{-1} \qquad;\qquad \oalign{\boldsymbol{L}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^{p} = \dot{\oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^{p}} {\oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^{p}}^{-1} = \displaystyle \sum_{s=1}^{M_{ss}} \dot{\gamma}^s \boldsymbol{m}^s \otimes \boldsymbol{n}^s\] where \(M_{ss}\) is the total number of slip systems, \(\dot{\gamma}^s\) denotes the plastic slip rate on the \(s\)th slip system, \(\boldsymbol{m}^s\) and \(\boldsymbol{n}^s\) are respectively the slip direction and slip plane normal of the \(s\)th slip system. Face-centered cubic (FCC) copper material is considered so that \(M_{ss}=12\) with \(\langle 110 \rangle\{111\}\) slip systems. The constitutive law is expressed here by the second Piola–Kirchoff stress tensor \[\oalign{\boldsymbol{\Pi}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \mathbb{C}:\oalign{\boldsymbol{E}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^{e}\] It follows that the Mandel stress tensor is given by \[\oalign{\boldsymbol{M}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = {\oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^{e}}^{T} \, \oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^{e} \, \oalign{\boldsymbol{\Pi}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\] For each slip system, plasticity is activated through a yield function \[f^s(\oalign{\boldsymbol{M}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}) = \rvert \tau^s\lvert - \tau^s_c\] where \(\tau^s\), is the resolved shear stress on the \(s\)th slip system and is related to the Mandel stress based on the relation \[\tau^s = \oalign{\boldsymbol{M}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} : \left(\boldsymbol{m}^s \otimes \boldsymbol{n}^s\right)\] and \(\tau_c^s\) denotes the critical resolved shear stress (CRSS) for the \(s\)th slip system, accounting for the yielding stress and a hardening. \(\tau_c^s\) evolves following a dislocation density-based hardening \[\tau^s_c = \tau^s_0 + \mu_h \, b \sqrt{ \displaystyle \sum_{u=1}^{M_{ss}} a^{su}\rho^u}\] with \(\tau^s_0\) the the thermal component of critical resolved shear stress of the \(s\)th slip system, \(\mu_h\) a hardening parameter, \(b\) the Burgers vector magnitude, \(a^{su}\) the interaction coefficient between dislocations on slip systems \(s\) and \(u\) and \(\rho^u\) the dislocation density on the slip system \(u\). In general, the same value of \(\tau^s_0\) is used for all the slip systems, i.e., \(\tau^s_0 = \tau_0\). The rate of the dislocation density on a given slip system \(s\) accounts for the multiplication and the annihilation through the formula \[\dot{\rho}^s = \displaystyle \frac{\lvert \dot{\gamma}^s\rvert}{b} \left( \displaystyle \frac{\sqrt{\displaystyle \sum_{u\neq s} \rho^u}}{\kappa} - yb \rho^s \right)\] with \(\kappa\) and \(y\) dimensionless coefficients which respectively characterize the dislocation mean free path and the dynamic recovery. A rate-dependent flow rule evolution is used \[\dot{\gamma}^s = \mathrm{sign}(\tau^s) \left(\max \left(\displaystyle \frac{f^s}{K_v},0\right)\right)^{n_{vp}}\] where the parameters \(K_v\) and \(n_{vp}\) govern viscosity. These constitutive equations have been implemented in the MFront code generator [39]. In this example, we choose the following values: for the self hardening \(a^{ss} = 0.1\), for the latent hardening \(a^{su} = 0.12, s\neq u\), for the viscosity parameters \(n_{vp}=30\), \(K_v=0.1\,\mathrm{MPa.s}^{1/n_{vp}}\) (note that the choice of the viscosity parameters leads a quasi-rate-independent behavior), for the critical resolved shear stress \(\tau_0=30\,\mathrm{MPa}\), for the hardening parameter \(\mu_h=9400\,\mathrm{MPa}\), for the dislocation density evolution, \(\kappa = 20\), \(b=2.54 \,\si{\angstrom}\), \(y=10\) and \(\rho_0^{tot} = 10^{12}\,\mathrm{m^{-2}}\). The elastic moduli are \(C_{11}=168\,\mathrm{GPa}\), \(C_{12}=121\,\mathrm{GPa}\) and \(C_{44}=75\,\mathrm{GPa}\).

[45] designed a torsion-bending experiment to study size effects on a \(\left< 111\right>\) oriented single crystal copper L-beam. Fig. 12 (a) presents an experimental setup example of a L-beam under a force-controlled actuator (applied force \(F\)). For the numerical simulations, the setup in Fig. 12 (b) is employed. Null Dirichlet BCs are applied on the left face \(S_{10}\) (normal to the axis \(\boldsymbol{e}_1\)) to mimic embedding. A Neumann BC is applied on the top face \(S_{21}\) (normal to the axis \(\boldsymbol{e}_2\)) with a non-null stress vector \(\boldsymbol{T^*}=-T^*_2 \boldsymbol{e}_2\), in a restricted half-circle region (i.e., an applied force \(F = T^*_2 \,A_{area}\) with the actuator area \(A_{area}=\pi a^2/8\), see Fig. 12 (b)). We applied a maximum level of stress of \(T^*_{2\,max}=60\,\mathrm{MPa}\) in \(3300\) loading steps. Furthermore, in order to have a unit-cell in a parallelepiped shape, the region outside the interested L-beam arm is filled with a void (i.e., an elastic material with null elastic moduli, see green part in Fig. 12 (b)).

Unlike in the previous section, \(\mathbb{C}_0\)-discrete Green operator can not be used for the present torsion-bending loading. In the context of such complex simulations, the optimal choice of pre-conditioner parameter(s), \(\mu_0\) (and \(\lambda_0\) for the \(\mathbb{B}_0\)-discrete Green operator) remains an open question. A first option is the classical Moulinec and Suquet [2] choice, i.e., \(\beta_0 = 1/2 \left(\max (\beta) + \min (\beta)\right)\) with \(\beta \in \{\lambda,\mu\}\). With this choice, the \(2\mu_0\)-discrete Green operator is outperformed by the \(\mathbb{B}_0\)-discrete Green operator. The same observation is made with regard to the choice \(\beta_0 = \max (\beta)\) with \(\beta \in \{\lambda,\mu\}\). In addition, the latter option results in a total of approximately \(37500\) iterations for the entire loading, as opposed to a total of \(45000\) iterations for the classical Moulinec and Suquet option, while employing the \(\mathbb{B}_0\)-discrete Green operator. For the purpose, we decide to continue all analysis with the \(\beta_0 = \max \left(\beta\right)\) and with the \(\mathbb{B}_0\)-discrete Green operator.

Fig. 13 (a) presents the macroscopic responses under the torsion-bending loading. Both finite difference schemes TETRA2 and HEX1 are used for comparison purpose. It turns out that TETRA2 scheme is once again more robust than HEX1 scheme. Actually the TETRA2 scheme allows to reach much larger plastic deformation level. This elastic-plastic response is qualitatively consistent with the experimental observations. Meanwhile, the HEX1 scheme fails to converge around the apparent yield stress.

Fig. 13 (b) displays the corresponding number of iterations required at each time step to satisfy the convergence criterion. After the very first loading step and for the most part of the elastic regime, the convergence iteration decreases with the loading so that a convergence on the equilibrium is reached with less than \(5\) iterations. The transition from elastic to plastic regime (micro-plasticity) leads to an increase in the required number of iterations for the convergence. It is shown that with HEX1 scheme a total of \(10^4\) iterations has not been enough to reach the convergence around the apparent yield stress.With the TETRA2, the convergence degrades around the yield stress but is reached with less than \(70\) iterations. In comparison with the HEX1 scheme, the TETRA2 scheme is demonstrated to produce a more stable evolution of the iterations from one loading step to another. Focusing on the results obtained with the TETRA2 scheme, within the plastic regime, a decrease of the convergence iteration is then observed (see Fig. Fig. 13 (b)).

The non-convergence of the HEX1 scheme at the apparent yield stress can be understood by visualizing the displacement field as shown in Fig. 14 (please note that the deformed shape is amplified for this particular Fig. 14 by a factor 10 only for visualization purpose). HEX1 leads to strong oscillations of the displacement field, preventing in the meantime a convergence as the loading increases. TETRA2 contributes to a smoother displacement field. Fig. 15 presents the deformed shape (this time without any amplification) of the L-beam in the present simulation (see Fig. 15 (a)) and in the experiment (see Fig. 15 (b)). In both cases, based on the geometrical design of the studied crystal, one can see that the arm in contact with the actuator is bent while the embedded arm is twisted, allowing to generate a strong strain gradient. A qualitatively good agreement with the experiment is observed.

In addition to the macroscopic behavior, it is worthwhile to examine the performance of the present MPI-based solver, for such simulations (finite transformation crystal plasticity with quasi-rate-independent). It is well known that these types of simulations are notably time-consuming (CPU time). Strong scalability tests are performed for the given discretization of \(n_1\times n_2 \times n_3 = 150\times50\times175\) voxels. The simulations are conducted on an in-house cluster that possesses \(10\) compute nodes with each node having the following characteristics: \(512\,\mathrm{Go}\) DDR5-RAM and \(2\) sockets with \(64\) cores per socket (AMD-EPYC9554) running at up to \(3.1\,\mathrm{GHz}\) frequency (i.e., a total of \(128\) cores per node). Fig. 16 (a) shows the raw computation time \(T_{\mathrm{simu}}\) as function of the number of CPU cores \(N_{\mathrm{CPUs}}\) (in practice, we increase the number of compute nodes from \(1\) to \(6\)). Considering the result with \(N_{\mathrm{CPUs}}^{\mathrm{ref}}=128\) CPU cores (i.e., a single node) as reference, the measured simulation time is compared to the ideal expected computation time, \(\dfrac{N_{\mathrm{CPUs}}^{\mathrm{ref}}\times T_{\mathrm{simu}}(N_{\mathrm{CPUs}}^{\mathrm{ref}})}{N_{\mathrm{CPUs}}}\). The measured time follows closely the ideal time. A significant decrease of the simulation time is observed. From \(128\) to \(640\) CPU cores, the simulation time drops from approximately from \(16\) hours to \(3.4\) hours while from \(640\) processors to \(768\) processors, the simulation time goes from \(3.4\) hours to \(3.2\) hours. Fig. 16 (b) presents these results in terms of the evaluation of the parallel efficiency, \(100\dfrac{N_{\mathrm{CPUs}}^{\mathrm{ref}}\times T_{\mathrm{simu}}(N_{\mathrm{CPUs}}^{\mathrm{ref}})}{N_{\mathrm{CPUs}} \times T_{\mathrm{simu}}(N_{\mathrm{CPUs}})}\). We obtain an excellent strong scaling maintaining approximately \(90\,\%\) efficiency up to \(640\) CPU cores. A drop in efficiency to \(79\,\%\) is observed at \(768\) processors. This could be the consequence of load imbalance inherent to the heterogeneous plastic strain distribution within the unit-cell. For the moment, the 2DECOMP&FFT library and AMITEX\(^\star\) consequently, only accounts for 2D pencils of fixed size during the simulation and load balance is not yet available. This point could be a future prospect for both 2DECOMP&FFT library and AMITEX\(^\star\) code.

4 Conclusion and recommendations↩︎

In this work, the FFT-based solver classically used to solve periodic mechanical problems [1], [2], is naturally extended to account for general (non-periodic and/or periodic) BCs, non-linear behaviors, small and finite transformations, within a parallel programming environment. The new developments are implemented within the AMITEX\(^\star\) code, which is a new version (in progress) of the open-source AMITEX_FTTP code [28]. These are based on new upgrades of the 2DECOMP&FFT library [30], [31], to account for discrete trigonometric transforms (DTTs).

The mechanical problem (small or finite transformation frameworks for Cauchy medium local theories) is reformulated with the fluctuation vector displacement as unknown. Leveraging between the link in the non-periodic BCs and the symmetry extensions of the fluctuation displacement, DTTs are introduced in a displacement-based accelerated fixed-point algorithm. The symmetry extensions of the concerned fields allow to build “virtual” extended periodic fields which discrete Fourier transform (DFT) can be related to the DTTs of non-extended fields. This helps to re-adapt the general aspects of the fixed-point algorithm with a minimal cost. Furthermore, finite difference schemes are used here to build the discrete Green operator. For this purpose, the recently introduced double tetrahedron [22], [27] and the classical hexahedral [24], [25] schemes are employed in the context of non-periodic BCs. Depending on the loading type and the small or finite transformation frameworks, three discrete Green operators, relying on different pre-conditioners (i.e., reference material), are discussed.

The present analysis draws upon the parallel programming concept to explore a range of numerical simulations inspired by realistic experimental non-trivial loading configurations. These simulations are employed to investigate the consequences of selecting a finite difference scheme and a discrete Green operator. We validated the presented approach against analytical solutions for elastic unit-cells loaded in bending (small strain) and in tension (finite transformation). For these examples, both the TETRA2 and HEX1 finite difference scheme results are in a very good agreement with the analytical solutions. Subsequent to this, isotropic plasticity and crystal plasticity material behaviors are employed within heterogeneous unit-cells. Based on the discussed cases (porous media, single-crystals), we show that for strong non-homogeneous loading in finite transformation framework, the TETRA2 scheme appears more robust than the HEX1 scheme. Concerning the discrete Green operator, its \(\mathbb{C}_0\)-based version converges faster than the two others, \(2\mu_0\)-based and \(\mathbb{B}_0\)-based versions, when the loading consists of “normal-mixed” BCs (compatible with the usage of \(\mathbb{C}_0\)-discrete Green operator). In other cases, the \(\mathbb{B}_0\)-discrete Green operator, is found to perform better. However, the optimization of the pre-conditioner parameters requires a deeper analysis, specially in the context of finite transformation. To conclude, the different applications demonstrated the versatility, robustness and parallel capabilities of the proposed FFT-based solver.

Going forward with the presented solver, it will be interesting to impose desired displacement and/or force values on interior nodes allowing to simulate for instance a compact tension specimen for which the loading is applied inside the unit-cell. Building on the various treated examples, it would be valuable, in a future work, to integrate the non-uniform grid procedures [8], [54], [55] and to enrich the finite difference schemes (e.g., with C3D20 FE type scheme) within the current generic FFT-based solver. Moreover, non-local models [47], [49], [56] could benefit from the versatile presented framework to allow for finite transformation scenarios and three-dimensional calculations, thereby enabling comparison with experimental high-resolution results. For multi-physic applications, it could also be interesting to discuss Robin-type BCs with the present framework.

Data availability↩︎

AMITEX\(^\star\) code used for the simulations, is available on request. It will be soon made open-source similarly to AMITEX_FTTP [28].

Acknowledgements↩︎

This research project was funded by the PTC-SN program (Programmes Transversaux de Competences-Simulation Numérique) of the French Alternative Energies and Atomic Energy Commission (CEA), France. Authors are grateful to Jean-Michel Scherer for providing the Mfront implementation of the crystal plasticity behavior law in finite transformation framework.

5 Symmetric gradient-conjugate \(\mathbb{C}_0\)-discrete Green operator↩︎

Following the detailed analysis for the expression of the discrete Green operator \({}_d\widehat{\mathrm{DGO}}(\widehat{\boldsymbol{p}})\) in Sec. 2.6.1 within finite transformation framework, we outline in the present appendix the changes within small strain framework. With small strain hypothesis, the pre-conditioner is considered as \[\mathbb{R}_0:{}_d\oalign{\boldsymbol{\nabla^s u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \mathbb{C}_0:{}_d\oalign{\boldsymbol{\nabla^s u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = 2\mu_0 \, {}_d\oalign{\boldsymbol{\nabla^s u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} + \lambda_0 \, \mathrm{tr}\left({}_d\oalign{\boldsymbol{\nabla^s u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\right) \boldsymbol{1}\] It implies that \[\widehat{\boldsymbol{p}} = \mu_0 \, \widehat{\boldsymbol{{}_d\nabla}\cdot\oalign{\boldsymbol{{}_d\nabla u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}} + (\mu_0 + \lambda_0) \, \widehat{\boldsymbol{{}_d\nabla}\cdot\oalign{\boldsymbol{{}_d\nabla^T u^f}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}}\] For a given finite difference scheme and for “normal-mixed” BCs, we show that the \({}_d\widehat{\mathrm{DGO}}(\widehat{\boldsymbol{p}})\) is given by

• HEX1 scheme: \[\widehat{\boldsymbol{u^{f}}} = - \displaystyle \frac{1}{\mu_0 \lVert {}_{H}\boldsymbol{\xi} \rVert^2} \left(\widehat{\boldsymbol{p}} - \frac{ \mu_0 + \lambda_0}{2\mu_0 + \lambda_0} \frac{\widehat{\boldsymbol{p}} \cdot {}_{H}\boldsymbol{\xi}}{ \lVert {}_{H}\boldsymbol{\xi} \rVert^2} {}_{H}\boldsymbol{\xi} \right) = {}_H\widehat{\mathrm{DGO}}(\widehat{\boldsymbol{p}})\]

• TETRA2 scheme: \[\widehat{\boldsymbol{u^{f}}} = \displaystyle - \frac{1}{\mu_0\,\lVert {}_{T_2}\boldsymbol{\xi} \rVert^2} \left[\widehat{\boldsymbol{p}} + \displaystyle \frac{ \mu_0 + \lambda_0}{2}\left(p_{\mathrm{int}1}\,\overline{{}_{T_2}\boldsymbol{\xi}} + p_{\mathrm{int}2}\,{}_{T_2}\boldsymbol{\xi} \right)\right] = {}_{T_1T_2}\widehat{\mathrm{DGO}}(\widehat{\boldsymbol{p}})\] with \(p_{\mathrm{int}1}\) and \(p_{\mathrm{int}2}\) intermediate fields defined as \[\left\{ \begin{array}{l} p_{\mathrm{int}1} = \displaystyle \frac{1}{D_s} \left[\displaystyle -\frac{3\mu_0+\lambda_0}{2}\,\lVert {}_{T_2}\boldsymbol{\xi} \rVert^2\,\left( \widehat{\boldsymbol{p}} \cdot {}_{T_2}\boldsymbol{\xi}\right) + \displaystyle \frac{\mu_0+\lambda_0}{2}\, \left({}_{T_2}\boldsymbol{\xi} \cdot {}_{T_2}\boldsymbol{\xi}\right) \left( \widehat{\boldsymbol{p}} \cdot \overline{{}_{T_2}\boldsymbol{\xi}}\right)\right] \\[1em] p_{\mathrm{int}2} = \displaystyle \frac{1}{D_s} \left[ \displaystyle \frac{\mu_0+\lambda_0}{2}\,\left(\overline{{}_{T_2}\boldsymbol{\xi}} \cdot \overline{{}_{T_2}\boldsymbol{\xi}}\right)\,\left( \widehat{\boldsymbol{p}} \cdot {}_{T_2}\boldsymbol{\xi}\right) - \displaystyle\frac{3\mu_0+\lambda_0}{2} \, \lVert {}_{T_2}\boldsymbol{\xi} \rVert^2 \left( \widehat{\boldsymbol{p}} \cdot \overline{{}_{T_2}\boldsymbol{\xi}}\right)\right] \\[1em] \qquad \text{where } D_s = \displaystyle \left(\frac{3\mu_0+\lambda_0}{2}\,\lVert {}_{T_2}\boldsymbol{\xi} \rVert^2\right)^2 - \left(\frac{\mu_0+\lambda_0}{2}\,\lvert {}_{T_2}\boldsymbol{\xi} \cdot {}_{T_2}\boldsymbol{\xi} \rvert\right)^2 \end{array} \right.\]

We recall that the modified frequencies \({}_{H}\boldsymbol{\xi}\) (real values), \({}_{T_1}\boldsymbol{\xi}\) and \({}_{T_2}\boldsymbol{\xi}\) (complex values) depend on the BCs.

6 Neo-Hooken beam under pure tension with Dirichlet BCs: numerical simulation versus (semi-)analytical finite transformation solution↩︎

This section proposes a validation of a displacement-controlled loading numerical simulation based on the presented FFT-based solver against an analytical solution. We consider an elasticity finite transformation framework to study a beam under a pure tension as sketched in Fig. 17 (a). In order to diversify the behavior law, a Neo-Hooken compressible elastic evolution is used in this example.

The strain energy density function \(\psi\) is given as \[\psi = \displaystyle \frac{\mu}{2} \left(\overline{I_1}-3\right) + \frac{K}{2}\left(J-1\right)^2\] where \(\mu\) and \(K\) are materials parameters respectively the shear and bulk moduli, \(\overline{I_1}\) and \(J\) are principal invariants defined as \[\left\{ \begin{array}{c} I_1 =\mathrm{tr}\left(\oalign{\boldsymbol{B}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\right) = \mathrm{tr}\left(\oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^T\right) \\[0.2em] I_2 = \displaystyle \frac{1}{2} \left(\mathrm{tr}\left(\oalign{\boldsymbol{B}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\right)^2 - \mathrm{tr}\left(\oalign{\boldsymbol{B}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}^2\right)\right) \\[0.2em] I_3 = \mathrm{det}\left(\oalign{\boldsymbol{B}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\right) = J^2 \end{array} \right. \Longrightarrow \left\{ \begin{array}{c} \overline{I_1} = \displaystyle J^{-2/3} I_1 \\[0.2em] \overline{I_2} = \displaystyle J^{-4/3} I_2 \end{array} \right.\] The Piola–Kirchhoff stress tensors are given by \[\oalign{\boldsymbol{\Pi}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \displaystyle \dfrac{\partial \psi}{\partial \oalign{\boldsymbol{E}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}} \quad; \quad \oalign{\boldsymbol{P}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} = \oalign{\boldsymbol{F}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth} \,\oalign{\boldsymbol{\Pi}\crcr\hidewidth\scriptscriptstyle\boldsymbol{\sim}\hidewidth}\] This behavior law is also implemented with MFront [39].

We demonstrate that the pure tension configuration leads to the analytical solution on the components \(P_{11}\) and \(P_{22}\) of the first Piola–Kirchhoff stress tensor in the form \[\left\{ \begin{array}{c} P_{11} = \displaystyle \frac{2}{3} \mu e_T^{2} \left( e\, e_T^{2} \right)^{-5/3} \left( e^{2} - e_T^{2} \right) + K \left( e\, e_T^{4} - e_T^{2} \right) \\[0.5em] P_{22} = \displaystyle \frac{1}{3} \mu \,e\, e_T \left( e\, e_T^{2} \right)^{-5/3} \left( e_T^{2} - e^{2} \right) + K e\, e_T \left( e_T^{2} - 1 \right) \end{array} \right.\] with \(e = 1 + \displaystyle \frac{u^*_1}{l}\) the longitudinal elongation and \(e_T\) the transversal elongation. The pure tension configuration results in the condition \(P_{22}=0\). Solving the non-trival equation \(P_{22} = 0\), helps to deduce the transversal elongation \(e_T\) from the elongation \(e\). In order to facilitate the calculation of \(e_T\), a simple Newton-Raphson method is employed in Python. Subsequently, a comparison can be drawn between the (semi-)analytical solution and the numerical simulations, depending on the finite difference scheme. Fig. 17 (b) presents the macroscopic responses. One can observe that the numerical simulations are in perfect accordance with the (semi-)analytical solution. Both schemes converge to the exact same non-linear responses.

References↩︎