# FFT-based disclination mechanics for high angle grain boundaries and comparisons with atomistic simulations

Abdallah Wazne, Joé Petrazoller, Julien Guénolé, Thiebaud Richeton, Stéphane Berbenni

► To cite this version:

Abdallah Wazne, Joé Petrazoller, Julien Guénolé, Thiebaud Richeton, Stéphane Berbenni. FFT-based disclination mechanics for high angle grain boundaries and comparisons with atomistic simulations. Mechanics Research Communications, 2025, 149, pp.104529. 10.1016/j.mechrescom.2025.104529. hal-05290382v2

HAL Id: hal-05290382

https://hal.univ-lorraine.fr/hal-05290382v2

Submitted on 3 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# FFT-based disclination mechanics for high angle grain boundaries and comparisons with atomistic simulations☆

Abdallah Wazne<sup>ID</sup>, Joé Petrazoller<sup>ID</sup>, Julien Guénolé<sup>ID</sup>, Thiebaud Richeton<sup>ID</sup>, Stéphane Berbenni<sup>ID*</sup>

Université de Lorraine, CNRS, Arts et Métiers Institute of Technology, LEM3, F-57000 Metz, France

## ARTICLE INFO

**Keywords:**

Grain boundaries

Disclinations

Anisotropic elasticity

FFT

Molecular statics

## ABSTRACT

The paper addresses the question of determining local stress fields using the numerically efficient Fast Fourier Transform (FFT) method as an application of the continuum defect theory in the presence of interfacial defects, and specifically High Angle Grain Boundaries (HAGBs). First, the Field Dislocation and Disclination Mechanics (FDDM) equations are reported highlighting the use of the Stokes-Helmholtz orthogonal decomposition for both elastic strain and curvature tensors which are involved for disclination-type defects. Then, following this decomposition, both incompatible and compatible fields can be solved. Second, the Green's function method for heterogeneous media is used to derive the stress polarization field, which integrates both dislocation and disclination densities. The compatible elastic strain tensor is numerically solved with the spectral method using the FFT (fast Fourier transform) algorithm together with finite difference (FD) schemes for computing first- and second-order spatial derivatives (FDDM-FFT). As applications of the method, the stress fields of two specific HAGBs (symmetric tilt GBs with [001] tilt axis), precisely $Σ29(520)[001]46.40°$ and $Σ149(1070)[001]20.02°$ are obtained from the Disclination Structural Unit Model (DSUM). They are calculated assuming both isotropic and anisotropic elasticity with the present FDDM-FFT numerical method and assuming periodic disclination density tensors. Quantitative comparisons are first performed for both HAGBs with analytical solutions obtained from specific combinations of disclination dipole walls in linear isotropic elasticity. Then, the effect of anisotropic elasticity is analyzed for both HAGBS considering two different FCC metals, namely Al and Ag. Lastly, some comparisons between the FDDM-FFT-based results with molecular statics (MS) simulations, using the virial stress method and an interpolation method based on Gaussian kernel are reported for both HAGBs applied to Al and Ag. It is shown that, despite their relative simplicity in describing HAGB defect cores, the FDDM-FFT results reproduce the major trends of MS-based results for both hydrostatic and shear stress components.

## 1. Introduction

Spectral methods that use the Fast Fourier Transform (FFT) algo-
rithm have seen a surge in interest in the micromechanics community
in recent decades [1,2]. Since the introduction of the FFT-based method
to solve micromechanical fields in composite materials as an alternative
to the Finite Element Method (FEM) by Moulinec and Suquet [3,4],
Lebensohn and co-workers have further developed FFT-based methods
for polycrystal viscoplasticity [5,6] and elastoviscoplasticity [7-9]. The
advantage of using the FFT formulation over the FEM was reported
in [10,11]. Subsequent numerical implementations of the FFT-based
method were developed for more and more complex local or non-local
constitutive behaviors. They include micromechanical problems with
eigenstrain fields [12-15], dilatational plasticity [16], strain-gradient

crystal plasticity [17–20], micropolar elasticity and plasticity [21,22],
phase field fracture [23,24], fatigue [25], mesoscale field dislocation
mechanics [26–30] and discrete dislocation dynamics (DDD) simula-
tions [31,32], providing the efficiency needed for the implementation
of these powerful and numerically demanding formulations. FFT-based
methods were adapted to compute elastic fields of continuum disloca-
tion mechanics [33–35] and generalized disclinations [34,36]. Interfa-
cial defects such as disconnection defects were modeled by developing
a field disclination mechanics theory based on FFT in [37,38].

Grain boundaries (GBs) are interfacial defects whose modeling in
terms of geometry (kinematics), elastic fields, and excess free energy
has motivated extensive research over the past century [39]. It is
well known that a GB has five independent macroscopic degrees of

<footnote>☆ This article is part of a Special issue entitled: 'Lebensohn SI_invited only' published in Mechanics Research Communications.</footnote>

<footnote>* Corresponding author.
E-mail address: stephane.berbenni@univ-lorraine.fr (S. Berbenni).</footnote>
freedom, i.e. three associated to the misorientation between the two
grains and two to define the grain boundary plane. In addition, four
microscopic degrees of freedom correspond to the translation vector
between both crystals and to the position of the GB along the normal
to the boundary plane. From the kinematical point of view, GBs are
referred to as tilt, twist, or mixed boundaries. Whatever its misorien-
tation angle, a tilt GB is defined by a rotation vector that is contained
within the boundary plane, while a twist GB is defined by a rotation
vector that is perpendicular to the GB plane. GBs are also classified as
Low-Angle Grain Boundaries (LAGBs) or High-angle grain boundaries
(HAGBs) depending on their misorientation angle that impact their
structure/excess energy relationship [40].

Based on the original work of Read and Shockley [41] and Read
[42], LAGBs are generally modeled by arrays of edge (resp. screw)
dislocations for tilt (resp. twist) GBs. Such models give satisfactory
results in terms of excess energy for GB with misorientations not
exceeding about 15°. Recently, Petrazoller et al. [43] calculated the
stress and strain fields of different LAGBs for applications to solute
segregation. These fields were evaluated by molecular statics (MS)
or by considering arrays of edge dislocations within the framework
of linear isotropic elasticity or anisotropic elasticity using the Stroh
formalism [44]. From atomistic simulations, the mechanical stresses
correspond to the symmetric virial stresses, where the per-atom stress
tensors [45] can be computed directly from interatomic forces derived
from interatomic potentials [46]. Elastic strains per-atom can also be
obtained by calculating the gradient of elastic displacements from an
ideal strain/stress-free crystal lattice [47].

For HAGBs, several models were developed in the literature. Among them, the geometrical model based on the coincident site lattice (CSL) model [48] was widely used to describe the grain boundary structure. It is based on the Σ value, which characterizes the ratio of the volumes of the coincidence unit cell over the original crystal lattice. Another well-known geometrical model is the Frank-Bilby model [49–52]. This model considers HAGBs as surface dislocations with a certain density resulting from the eigendistortion due to GB dislocations [51,53]. The Structural Unit Model (SUM) [54–57] is an extension of the Read and Shockley’s model. In this model, the GB structure is described as a combination of building blocks called structural units (SUs), each having a well-defined atomic structure, rather than as an array of dislocations. More recently, it has been reported that GB metastability is important to characterize GB structure and energy in relation to statistical mechanics of GB states [58]. Then, for a given GB type, there is a broad spectrum of energies [58], which can be explored to generalize the SUM [59].

An alternative approach to modeling HAGBs involves identifying SUs from atomistic simulations and bridging the gap between the discrete and continuum representation of GBs using the Disclination Structural Unit Model (DSUM). Atomistic simulations such as MS or molecular dynamics (MD) are indeed well suited to identify internal structure, i.e. intrinsic defects within LAGBs or HAGBs [52,60–62], and interface-driven mechanisms at the nanoscale [63–69]. These approaches provide unique access to the discrete static and dynamic atomic scale configurations of the defects, with assumptions mostly limited to atomic interactions.

On the other hand, the DSUM model has been widely applied to
determine the elastic fields and the energy of GBs [70–72]. Since
the original work of Li [73], disclinations have been introduced to
describe grain and subgrain boundaries, twin junctions, and deforma-
tion defects in nanocrystalline materials, which were often represented
as arrays of disclination dipoles [69,74–79]. They were also used to
model disconnection defects [38,80,81]. The theoretical formulation
of dislocations and disclinations in an elastic medium was initially
developed by Volterra [82], where these defects were described as
displacement and rotational discontinuities, respectively. When disclia-
nations coexist with dislocations in crystalline materials, both the
elastic strain and the curvature tensors become incompatible, resulting

in multi-valued displacement and rotational fields [53,83–88]. The
concept of a disclination density tensor was introduced in [83,85,89],
and subsequent developments established a theoretical framework for
disclination dynamics [90]. Recent advances led to the Field Disloca-
tion and Disclination Mechanics (FDDM) theory, which extends field
dislocation mechanics (FDM) to capture non-local interactions arising
from both defect types [79,91–93]. A key distinction of FDDM com-
pared to more classic eigenstrain-based micromechanics [12,14,15,94]
is the Stokes-Helmholtz decomposition of elastic fields into compatible
and incompatible components [95]. This decomposition provides a
rigorous basis for quantifying elastic incompatibilities induced by bulk
and interfacial defects. For example, atomistic simulations have been
recently used in combination with FDM to explore the atomic scale
origin of strain-gradient approaches modeling grain boundaries [96].

The main objective of the present paper is to calculate the elastic
fields (essentially stresses) of two symmetric tilt HAGBs with [001]
tilt axis, precisely Σ29(520)[001]46.40° and Σ149(1070)[001]20.02°.
This choice of both HAGBs was made because they are well-defined
stable [001] symmetric tilt GBs constituted of simple combinations of
well-known SUs. In comparison with our previous work [43], where
only LAGBs were studied, the present study investigates HAGBs with
misorientation being either a little higher than the classic limit of
LAGBs (estimated to approximately 15°) for Σ149, or much higher than
this limit for Σ29. After having identified their SUs, both HAGBs can
be easily described by the DSUM with two different combinations of
wedge disclination dipoles (line defects). They are first built up from
atomistic simulations and then with FFT-based field dislocation and
disclination mechanics (FDDM) using their DSUM representation. The
disclination density tensor distribution in the DSUM for wedge disclina-
tion dipoles is set up with calibrated Gaussian distributions resolved on
the FFT grid that are validated in comparison with analytical solutions
in linear isotropic elasticity. Then, the role of elastic anisotropy is
highlighted by comparing the results for both HAGBs of two cubic
materials with different Zener ratio: Al and Ag. Finally, the predictions
of the FFT-based DSUM model in anisotropic elasticity are compared
to those obtained from atomistic results.

The paper is organized as follows. First, the notations are explained
at the end the present introduction, providing the necessary definitions
and symbols used throughout the study. Section 2 details the Fast
Fourier Transform (FFT)-based method for solving the elasto-static field
equations of field disclination mechanics (FDDM). Section 3 introduces
the atomistic construct of both studied HAGBs, establishing a con-
nection between atomistic simulations and the DSUM. In Section 4,
various study cases are examined, including the stress fields of HAGBs
using FFT versus analytical solutions in isotropic elasticity, the FFT
predictions of stress fields for the same HAGBs in Al and Ag, com-
paring anisotropic and isotropic elasticity, and finally a comparison of
stress field predictions between FFT-based anisotropic modeling and
atomistic simulations. Finally, Section 5 outlines the conclusions and
perspectives of this work.

Notations

Bold symbols are used to represent tensors or vectors. When an
index notation is used, the Einstein's summation convention on dummy
indices is applied. The symmetric part of a second-order tensor **A** is
denoted as **A**<sup>*sym*</sup>. When expressed in components, a symmetric tensor
is written as **A**<sub>*ij*</sub>, where **A**<sub>*ij*</sub> = **A**<sub>*ji*</sub>. The transpose of a second-order
tensor **A** is represented as **A**<sup>*t*</sup>, with components **A**<sup>*t*</sup><sub>*ij*</sub> = **A**<sub>*ji*</sub> in indicial
notation. The dot product of two tensors, **A** · **B**, produces a tensor with
components **A**<sub>*ik*</sub>**B**<sub>*kj*</sub>, while their tensor product **A**⊗**B** yields components
**A**<sub>*ij*</sub>**B**<sub>*kl*</sub>. Similarly, the dot product of a tensor **A** with a vector **V** results
in a vector **A** · **V**, whose components are **A**<sub>*ij*</sub>**V**<sub>*j*</sub>. The trace inner product
of two second-order tensors is denoted as **A** : **B** = **A**<sub>*ij*</sub>**B**<sub>*ji*</sub>, or, in the case
of a higher-order tensor and a second-order tensor, **A** : **B** = **A**<sub>*ijkl*</sub>**B**<sub>*kl*</sub>.
<figure><img src="image_4.png" /><figcaption>Fig. 1. Stable atomic structures of Al Σ29(520)[001]46.40° (a), and Al Σ149(1070)[001]20.02° HAGBs, with the identification of their minor and major structural units (SUs). Atoms belonging to the same plane are represented with the same color. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 2. Disclination Structural Unit Model (DSUM) representation of (a) Σ29(520)[001]46.40° and (b) Σ149(1070)[001]20.02° symmetric tilt grain boundaries.</figcaption></figure>

In rectangular Cartesian coordinates, the operations involving the cross product “×” of a second-order tensor A and a vector V, the gradient (grad) of a vector V, as well as the divergence (div) and curl (curl) of A, are expressed as follows:

$$ (\mathbf{A} \times \mathbf{V})_{ij} = e_{jkl} A_{ik} V_l $$

$$ (\mathbf{grad} \mathbf{V})_{ij} = V_{j,i} $$

$$ (\mathbf{div} \mathbf{A})_i = A_{ji,j} $$

$$ (\mathbf{curl} \mathbf{A})_{il} = e_{imk} A_{kl,m} $$

Here, $e_{pmk}$ represents the components of the third-order Levi-Civita tensor, and spatial derivatives with respect to Cartesian coordinates are indicated by a comma followed by the respective component index.

## 2. Fast fourier transform-based defect mechanics

### 2.1. Field dislocation and disclination mechanics (FDDM)

Following the continuum theory of defects introduced by Kröner [97,98] for dislocations and by De Wit [85,99] for more general defects

including disclinations, let us introduce the elastic curvature tensor $\kappa^e$ and the elastic strain tensor $\epsilon^e$ in a body V that contains defects of both disclination and dislocation types. We first define a closed circuit C (generally called *Burgers circuit*) delimiting a surface S of unit normal n so that disclination and dislocation lines cross the surface S. In the case where the resulting Frank and Burgers vectors of all threading disclination/dislocation lines are non-zero, there is incompatibility, in the sense that the rotation and displacement fields become multi-valued. A measure of this incompatibility for a simply connected body is precisely the net Frank vector $\omega$ and the net general Burgers vector b respectively. Both quantities are defined along the Burgers circuit C using Weingarten's theorem [85]:

$$ \omega = \int_C \kappa^e \cdot dl \qquad (1) $$

$$ b = \int_C (\epsilon^e - \kappa^e \times \mathbf{x}) \cdot dl $$
<figure><img src="image_6.png" /><figcaption>Fig. 3. Interpolation of $\sigma_{kk}$ from MS and FDDM-FFT calculations of the Ag $\Sigma29(520)[001]$46.40° boundary, using a Gaussian law with $\bar{\sigma} = 1.2$ Å for atomistic results, and $\bar{\sigma} = 0.15$ Å for FFT results. The interpolation cutoff is set to $3\bar{\sigma}$ (see Eqs. (25)).</figcaption></figure>

where a FS/RH convention is used. Using Stokes' theorem, Eq. (1) yields:

$$ 
\begin{aligned}
\omega &= \int_S \theta \cdot \mathbf{n} dS \\
\mathbf{b} &= \int_S (\alpha - \theta \times \mathbf{x}) \cdot \mathbf{n} dS
\end{aligned}
\qquad (2) 
$$

where the disclination density tensor $\theta$ and the dislocation density tensor $\alpha$ are introduced. Using Cartesian coordinates ($e_1, e_2, e_3$), their components $\theta_{mi}$ and $\alpha_{mi}$ represent the net Frank vector and Burgers vector, respectively, in the direction $e_i$ per unit area of the surface S, with the unit line vector oriented along $e_m$. The linear theory of static disclinations developed by De Wit [85] yields the following definitions of $\theta_{mi}$ and $\alpha_{mi}$:

$$ 
\begin{aligned}
\theta_{mi} &= e_{mjk} \kappa_{kl,j}^{e} \\
\alpha_{mi} &= e_{mjk} (\varepsilon_{(kl),j}^{e} + e_{ipk} \kappa_{ip}^{e})
\end{aligned}
\qquad (3) 
$$

Then, from Eqs. (3), conservation equations hold:

$$ 
\begin{aligned}
\theta_{ji,j} &= 0 \\
\alpha_{ji,j} + e_{ikm} \theta_{km} &= 0
\end{aligned}
\qquad (4) 
$$

From a physical standpoint, the first equation implies that disclinations cannot terminate within the body, while the second equation indicates that dislocations can only end on disclinations, as discussed in [86]. Going beyond De Wit [85]'s equations and following the previous work of Berbenni and Taupin [37], the elastic strain tensor $\varepsilon^e$ can be represented through the Stokes-Helmholtz orthogonal decomposition [95,100,101], where unique tensor fields $\chi$ and $w$ do exist such that:

$$ \varepsilon^e = \varepsilon^{e,\perp} + \varepsilon^{e,\|} = (\mathrm{curl} \, \chi)^{\mathrm{sym}} + (\mathrm{grad} \, w)^{\mathrm{sym}}, \qquad (5) $$

The orthogonality condition for the decomposition is given by: $\int_V \mathrm{curl} \chi : \mathrm{grad} \, w dV = 0$ where $\varepsilon^{e,\perp} = (\mathrm{curl} \, \chi)^{\mathrm{sym}}$ and $\varepsilon^{e,\|} = (\mathrm{grad} \, w)^{\mathrm{sym}}$

represent the incompatible and compatible components of $\varepsilon^e$, respectively. Similarly, the Stokes-Helmholtz orthogonal decomposition of the elastic curvature tensor $\kappa^e$ yields:

$$ \kappa^e = \kappa^{e,\perp} + \kappa^{e,\|} = \mathrm{curl} \, \psi + \mathrm{grad} \, z, \qquad (6) $$

with the orthogonality condition $\int_V \mathrm{curl} \, \psi : \mathrm{grad} \, z dV = 0$ where $\kappa^{e,\perp} = \mathrm{curl} \, \psi$ and $\kappa^{e,\|} = \mathrm{grad} \, z$ denote the incompatible and compatible components of $\kappa^e$, respectively. From $\mathrm{div} \, \mathrm{curl} \, \chi = 0$ and $\mathrm{div} \, \mathrm{curl} \, \psi = 0$, it comes:

$$ \mathrm{div} \, \kappa^{e,\perp} = 0 \qquad (7) $$

$$ \mathrm{div} \, \varepsilon^{e,\perp} = 0 \qquad (7) $$

By using the Stokes-Helmholtz decompositions (Eqs. (5) and (6)), and noting that $\mathrm{curl} \, \mathrm{grad} \, z = 0$ and $\mathrm{curl} \, \mathrm{grad} \, w = 0$, Eqs. (3) yield:

$$ 
\begin{aligned}
\theta_{mi} &= e_{mjk} \kappa_{kl,j}^{e,\perp} \\
\alpha_{mi} &= e_{mjk} (\varepsilon_{(kl),j}^{e,\perp} + e_{ipk} \kappa_{ip}^{e,\perp})
\end{aligned}
\qquad (8) 
$$

and:

$$ 
\begin{aligned}
e_{mjk} \kappa_{kl,j}^{e,\|} &= 0 \\
e_{mjk} (\varepsilon_{(kl),j}^{e,\|} + e_{ipk} \kappa_{ip}^{e,\|}) &= 0
\end{aligned}
\qquad (9) 
$$

Additionally, the compatible fields $\varepsilon^{e,\|}$ and $\kappa^{e,\|}$ yield:

$$ 
\begin{aligned}
\varepsilon_{(ij)}^{e,\|} &= \frac{1}{2} (w_{i,j} + w_{j,i}) \\
\kappa_{ij}^{e,\|} &= z_{j,i} = \frac{1}{2} e_{jkl} w_{l,ki}
\end{aligned}
\qquad (10) 
$$

Using $\mathrm{curl} \, \mathrm{curl} = \mathrm{grad} \, \mathrm{div} - \Delta$, where $\Delta$ is the Laplacian operator, then applying Eqs. (7) together with the identity: $e_{js}e_{jk} = \delta_{sk}\delta_{li} - \delta_{sl}\delta_{ik}$, $\kappa^{e,\perp}$ and $\varepsilon^{e,\perp}$ satisfy the following Poisson equations:

$$ 
\begin{aligned}
\kappa_{j,i,k,k}^{e,\perp} &= -e_{jkl} \theta_{li,k} \\
\varepsilon_{(lr),kk}^{e,\perp} &= -e_{rql} (\alpha_{il,q} - \delta_{il} \kappa_{pp,q}^{e,\perp} + \kappa_{il,q}^{t,e,\perp})_{(lr)}
\end{aligned}
\qquad (11) 
$$
<figure><img src="image_7.png" /><figcaption>Fig. 4. Comparison of the contour plots of the trace of the stress tensor $\sigma_{kk}$ and shear stress $\sigma_{xy}$ normalized by bulk modulus $K$ and $2\mu$ respectively, for a $\Sigma29(520)[001]46.40^\circ$ tilt grain boundary modeled using DSUM, obtained via FDDM-FFT and analytical methods [75]. Results were calculated with 1024 $\times$ 1024 pixels.</figcaption></figure>

where $\kappa_{ij}^{\tau e, \perp} = \kappa_{ij}^{e, \perp}$, and $\kappa_{pp}^{e, \perp}$ is the trace of $\kappa^{e, \perp}$. In the absence of body forces and inertia effects, a uniform stress $\bar{\sigma}$ is applied as the traction boundary condition. This corresponds to the spatial average of the stress field $\langle\sigma\rangle$ over the unit cell $V$, as per the spatial averaging theorem. The balance of linear momentum yields:

$$ \mathrm{div} \, \sigma = 0 \text{ in } V. \qquad (12) $$

We consider the generalized Hooke's law of a linear elastic anisotropic heterogeneous medium:

$$ \sigma = C : \varepsilon^e, \qquad (13) $$

where $C$ is the fourth-order tensor representing spatially heterogeneous elastic moduli, which exhibits classic symmetries for Cartesian components: $C_{ijkl} = C_{jikl} = C_{ijkl} = C_{klij}$. Using Eq. (5) in conjunction with Eq. (13), the balance of linear momentum (Eq. (12)) can be rewritten as a heterogeneous Navier-type equation in $V$, incorporating both compatible and incompatible fields:

$$ \mathrm{div} (C : \varepsilon^{e, \parallel}) + \mathbf{f}^{\perp} = 0, \qquad (14) $$

The fictitious body force density due to elastic incompatibilities (dislocations and dislocations) is $\mathbf{f}^{\perp} = \mathrm{div}(C : \varepsilon^{e, \perp})$. Considering a homogeneous reference medium characterized by linear elastic stiffness tensor

$C^0$, where $C(x) = C^0 + \delta C(x)$ (with $x$ denoting a position vector), Eq. (14) can be rewritten in indicial form:

$$ C_{ijkl}^{0} w_{l,kj} + \tau_{ij,j} = 0, \qquad (15) $$

where $\tau$ is the stress polarization tensor such that:

$$ \sigma = C^0 : \varepsilon^{e, \parallel} + \tau \qquad (16) $$

with:

$$ \tau = \sigma^{\perp} + \delta C : \varepsilon^{e, \parallel} \qquad (17) $$

In Eq. (17), the first term, $\sigma^{\perp} = C : \varepsilon^{e, \perp}$, represents the stress resulting from incompatible elastic strains, which are determined by solving the incompatibility problem alone (as given in Eqs. (11)). Since $\tau$ includes the unknown compatible elastic strain, $\varepsilon^{e, \parallel}$, Eq. (15) is solved by applying the Green's function method to end with the Lippmann-Schwinger's type equation for the unknown compatible elastic strain, $\varepsilon^{e, \parallel}$:

$$ \varepsilon^{e, \parallel}(\mathbf{x}) = \langle \varepsilon^{e, \parallel} \rangle - (\Gamma^0 \star \tau)(\mathbf{x}), \qquad (18) $$

Here, $\star$ denotes the spatial convolution product, and $\Gamma^0$ represents the strain-based Green's operator associated with the homogeneous reference elastic medium $C^0$. The term $\langle \varepsilon^{e, \parallel} \rangle$ indicates the spatial
<figure><img src="image_8.png" /><figcaption>Fig. 5. Comparison of the contour plots of the trace of the stress tensor $\sigma_{kk}$ and the shear stress $\sigma_{xy}$ normalized by bulk modulus $K$ and $2\mu$ respectively, for a $\Sigma149(1070)[001]20.02^\circ$ tilt grain boundary modeled using DSUM, obtained via FDDM-FFT and analytical methods [75]. Results were calculated with 1024 x 1024 pixels.</figcaption></figure>

average of $\epsilon^{e,\|}$ over the unit cell $V$ which is determined by averaging Eq. (16):

$$ \langle \epsilon^{e,\|} \rangle = C^{0-1} : (\bar{\sigma} - \langle \tau \rangle) \qquad (19) $$

where $\langle \tau \rangle$ is the average stress polarization over the unit cell $V$. Eqs. (11) together with Eqs. (17)–(19) are now numerically solved with a fast Fourier transform (FFT) method.

## 2.2. Fast fourier transform (FFT)-based numerical implementation: FDDM-FFT

The incompatible fields $\kappa^{e,\perp}$ and $\epsilon^{e,\perp}$ along with the compatible field $\epsilon^{e,\|}$, are determined in Fourier space. The corresponding elastic fields are subsequently reconstructed in real space through the application of the inverse Fourier Transform. In Fourier space, let $\xi$ represent the Fourier vector, where its magnitude is defined as $\xi = \sqrt{\xi^2 \cdot \xi}$, and its Cartesian components are denoted as $\xi_i$. The imaginary unit is expressed as $i = \sqrt{-1}$. Furthermore, let $\hat{\theta}(\xi)$, $\hat{\alpha}(\xi)$, $\hat{\kappa}^{e,\perp}(\xi)$, and $\hat{\epsilon}^{e,\perp}(\xi)$ be the Fourier transforms of the corresponding real-space fields $\theta(x)$, $\alpha(x)$, $\kappa^{e,\perp}(x)$, and $\epsilon^{e,\perp}(x)$. Both Poisson-type equations (Eqs. (11)) are

solved in Fourier space using the differentiation theorem:

$$ \begin{aligned} \hat{\kappa}_{ji}^{e,\perp}(\xi) &= i \frac{\xi_k}{\xi^2} e_{jkl} \hat{\theta}_{ll}(\xi) \quad \forall \xi \neq 0, \\ \hat{\kappa}_{ji}^{e,\perp}(\mathbf{0}) &= \mathbf{0}, \\ \hat{\epsilon}_{(lr)}^{e,\perp}(\xi) &= i \frac{\xi_q}{\xi^2} e_{rqi} \left( \hat{\alpha}_{il}(\xi) - \delta_{il} \hat{\kappa}_{pp}^{e,\perp}(\xi) + \hat{\kappa}_{il}^{r e,\perp}(\xi) \right)_{(lr)} \quad \forall \xi \neq 0, \\ \hat{\epsilon}_{(lr)}^{e,\perp}(\mathbf{0}) &= \mathbf{0}, \end{aligned} \qquad (20) $$

Let $\hat{\epsilon}^{e,\|}(\xi)$ and $\hat{\Gamma}^0(\xi)$ be the continuous Fourier transforms of $\epsilon^{e,\|}(x)$ and $\Gamma^0(x)$, respectively. By applying the Fourier transform to the integral Lippmann–Schwinger–Dyson type equation (Eq. (18)), we obtain:

$$ \begin{aligned} \hat{\epsilon}^{e,\|}(\xi) &= -\hat{\Gamma}^0(\xi) : \hat{\tau}(\xi) \quad \forall \xi \neq 0, \\ \hat{\epsilon}^{e,\|}(\mathbf{0}) &= \langle \epsilon^{e,\|} \rangle. \end{aligned} \qquad (21) $$

The Fourier transform of the modified Green operator $\hat{\Gamma}^0$ for anisotropic materials can be determined in Fourier space, see e.g. [53]. In the case of isotropic elasticity, the components of $\hat{\Gamma}^0(\xi)$ were given by Moulinec and Suquet [4].

The numerical method used to compute the balanced elastic fields associated with dislocation and dislocation density distributions in periodic media consists of two main steps: (i) the initialization of the
<figure><img src="image_9.png" /><figcaption>Fig. 6. Comparison of the contour plots of the trace of the stress tensor $\sigma_{kk}$ and shear stress $\sigma_{xy}$ for a $\Sigma29$ tilt grain boundary modeled using DSUM, obtained via FDDM-FFT for isotropic and anisotropic Al. Results were calculated with 1024 $\times$ 1024 pixels.</figcaption></figure>

iterative process, which involves setting $\epsilon^{e,\|}$ for a macroscopic applied stress $\bar{\sigma}$ after calculating $\hat{\kappa}^{e,\perp}(\xi)$ and $\hat{\epsilon}^{e,\perp}(\xi)$ using Eqs. (20), and (ii) the global iterative procedure, which follows a fixed-point algorithm (i.e. basic algorithm) to solve Eq. (21), where $\epsilon^{e,\|}$ is computed once convergence is achieved. Technical details about the FFT algorithm can be found in Appendix A.

## 2.3. Implementation of disclination dipoles on the FFT grid

For the present FFT calculations for HAGBs, the disclination density tensor $\theta$ is provided as the only source of defects in the form of disclination dipoles. For a straight wedge disclination, both the Frank and line vectors are along the z-axis. Therefore, the only non-zero prescribed disclination density tensor component is $\theta_{33}$. From Eqs. (2) and for a positive disclination, the latter is related to the Frank's vector magnitude $\omega$ as follows:

$$ \omega = \omega_3 = \int_S \theta_{33} dS \qquad (22) $$

For the numerical implementation, $\theta_{33}$ will be given by a two-dimensional radial normal Gaussian distribution for disclination dipoles (with $\pm$ signs):

$$ \theta_{33} = \pm \frac{\omega}{2\pi\sigma^2} \exp\left(-\frac{r^2}{2\sigma^2}\right) \quad \text{if } r < d \qquad (23) $$

$$ \theta_{33} = 0 \qquad \text{if } r \ge d $$

where $d$ is the FFT computational grid size, $r = \sqrt{x^2 + y^2}$ is the polar radius and $\sigma$ is the standard deviation of the Gaussian function. Specific sequences of disclination dipoles will be set up in the next section by modeling two specific HAGBs with the DSUM.

# 3. High Angle Grain boundary (HAGB): from atomistics to disclination structural unit model

## 3.1. HAGB construct from atomistics

To enable comparisons with the FFT method discussed in the previous section, reference configurations were generated using molecular statics (MS). This atomistic simulation technique relies on the minimization of an interatomic potential to determine the lowest-energy configuration of an atomic system. It uses optimization schemes for atomic positions that nullify net interatomic forces and reach a metastable location in the potential energy landscape [102]. Intrinsically, such approach is static: it cannot capture the physics of particle motion and does not consider temperature.

Two symmetric tilt HAGBs with the tilt axis oriented along the [001] crystallographic direction were considered: $\Sigma29(520)[001]$46.40° and $\Sigma149(1070)[001]$20.02°, within two different FCC materials (Al and Ag). These HAGBs were first constructed based on the CSL theory using the Python package GB_code [103]. Subsequently, using the embedded atom model (EAM) interatomic potentials of Liu et al. [104]
<figure><img src="image_10.png" /><figcaption>**Fig. 7.** Comparison of the contour plots of the trace of the stress tensor $\sigma_{kk}$ and shear stress $\sigma_{xy}$ for a $\Sigma29$ tilt grain boundary modeled using DSUM, obtained via FDDM-FFT for isotropic and anisotropic Ag. Results were calculated with 1024 $\times$ 1024 pixels.</figcaption></figure>

<table><caption>**Table 1**

Values of the lattice parameter *a*, cubic elastic constants and Zener ratio *A*, as computed with the two EAM potentials used in this work.</caption><thead><tr><th>System</th><th>Potential</th><th>a (Å)</th><th>C<sub>11</sub> (GPa)</th><th>C<sub>12</sub> (GPa)</th><th>C<sub>44</sub> (GPa)</th><th>A</th></tr></thead><tbody><tr><td>Al</td><td>Liu et al. [104]</td><td>4.05</td><td>118.9</td><td>62.9</td><td>32.9</td><td>1.22</td></tr><tr><td>Ag</td><td>Pan et al. [105]</td><td>4.16</td><td>123.4</td><td>94.1</td><td>47.1</td><td>3.21</td></tr></tbody></table>

and Pan et al. [105], for Al and Ag respectively (see Table 1), in-plane translations within the grain boundary plane were performed to explore microscopic degrees of freedom and identify the configuration with minimum energy. Energy minimization for each system was performed using the LAMMPS software (Version 2022-06-23 Update 1), employing the FIRE algorithm [106,107] with a force convergence criterion of $10^{-9}$ eV/Å. After identifying the lowest-energy configuration, the simulation box was relaxed to minimize residual stresses using the conjugate gradient algorithm with the same force criterion. A stable structure was then obtained, on which it was possible to identify minor and major structural units (SUs). These ones are of types C and B for the $\Sigma29$ HAGB and of types B and A for the $\Sigma149$ HAGB, see Figs. 1(a) and (b), respectively. Per-atom stress tensors [45] were computed directly from interatomic forces in units of pressure $\times$ volume, then converted into stresses using each atom's associated Voronoi volume [46].

### 3.2. HAGB construct from the disclination structural unit model (DSUM)

In parallel to the construct from atomistics, both HAGBs $\Sigma29(520)$ [001]46.40° and $\Sigma149(1070)$[001]20.02° were built up using the disclination structural unit model (DSUM). The DSUM was initially proposed by Li in 1972 [73] and later enhanced by Gertsman et al. in 1989 [70]. It creates a grain boundary with a specific misorientation angle $\theta$ by breaking it down into a contiguous and alternating sequence of M majority ($m$) and N minority ($n$) structural units (SUs). The $m$ and $n$ units have misorientation angles $\theta_m$ and $\theta_n$, respectively, with $\theta_m < \theta_n$. In the DSUM, the GBs are represented as a wall of wedge disclination dipoles bounding the minority SUs, see Figs. 2(a) and (b). For each dipole, the disclination strength or Frank's vector magnitude is defined as [71]: $\omega = \theta_n - \theta_m$ (for positive pole). The period of the GB is given by $H = ML_m + NL_n$, where $L_n$ represents the characteristic length of the minority SU with disclination dipole arm length equal to $0.5L_n$, and $L_m$ the characteristic length of the majority SU. These dimensions can be expressed through the following geometric relations [72]:

$$L_m = D_m \cos\left(\frac{\theta - \theta_m}{2}\right) \qquad L_n = D_n \cos\left(\frac{\theta_n - \theta}{2}\right) \qquad (24)$$

where $\theta$ is the average GB misorientation. The lengths $D_m$ and $D_n$ depend on the type of majority and minority SUs coming from the structural decomposition of the period, as expressed with respect to the lattice parameter *a* for A, B, C and D SUs: $D_A = a\frac{\sqrt{2}}{2}$, $D_B = a\frac{\sqrt{5}}{2}$,
<figure><img src="image_11.png" /><figcaption>**Fig. 8.** Comparison of the contour plots of the trace of the stress tensor $\sigma_{kk}$ and shear stress $\sigma_{xy}$ for a $\Sigma$149 tilt grain boundary modeled using DSUM, obtained via FDDM-FFT for isotropic and anisotropic Al. Results were calculated with 1024 $\times$ 1024 pixels.</figcaption></figure>

$D_C = a\sqrt{\frac{10}{2}}$ and $D_D = \frac{a}{2}$. The misorientation angles of these SUs are given by: $\theta_A = 0^{\circ}$, $\theta_B = 36.87^{\circ}$, $\theta_C = 53.13^{\circ}$ and $\theta_D = 90^{\circ}$.

As seen in Fig. 1(a), the $\Sigma$29 HAGB has a structural decomposition $|BC.BC|$ as reported in Table 2, where B is the majority unit and C is the minority unit. The periodicity is $H = L_m + L_n$ (see Fig. 2(a)). In contrast, the structural decomposition of the $\Sigma$149 boundary is $|AABABAB.AABABAB|$ (Fig. 1(b)) as reported in Table 2, where A is the majority unit and B is the minority unit. Its periodicity is given by $H = 4L_m + 3L_n$ (see Fig. 2(b)). The analytical stress field in isotropic elasticity of the $\Sigma$29 HAGB was modeled as a single wall of disclination dipoles with arm length $0.5L_n$ and strength $\pm\omega = \pm(\theta_n - \theta_m) = \pm 16.26^{\circ}$ (see Fig. 2(a)). For the $\Sigma$149 HAGB, it was modeled as the superposition of three infinite walls of disclination dipoles with arm length $0.5L_n$ and strength $\pm\omega = \pm(\theta_n - \theta_m) = \pm 36.87^{\circ}$. The first wall is centered at $y = 0$, the second at $y = L_m + L_n$, and the third at $y = 2(L_m + L_n)$ (see Fig. 2(b)). The analytical formula for the stress field of a single wall of disclination dipoles with period $H$, arm length $0.5L_n$, and strength $\omega$, is recalled in Appendix B.

For FFT calculations (see Section 2), the sequences of disclination dipoles were located at distances $\pm 0.5L_n$ from the center of each dipole in both HAGBs. The computational grid consisted of 1024 $\times$ 1024 pixels, with $d = H/128$ for the $\Sigma$149 HAGB and $d = H/64$ for the $\Sigma$29 HAGB.

<caption>**Table 2**</caption>

Structural unit decomposition of the two studied high angle symmetric tilt GBs about the $[001]$ axis [55,69]. $a$ is the lattice parameter. Bars denote one period of the boundary and the dot signifies that equivalent atoms in each half period are relatively displaced by $1/2[001]$ along the tilt axis [69].

<table><thead><tr><th>θ [°]</th><th>GB plane</th><th>Σ</th><th>SU period decomposition</th><th>Period vector</th></tr></thead><tbody><tr><td>46.40</td><td>(5 2 0)</td><td>29</td><td>|BC.BC|</td><td>a[2 5 0]</td></tr><tr><td>20.02</td><td>(10 7 0)</td><td>149</td><td>|AABABAB.AABABAB|</td><td>a[7 10 0]</td></tr></tbody></table>

### 3.3. Comparison method between FFT and atomistic calculations

The lattice parameters and the elastic constants considered in the FFT calculations are those of the interatomic potentials (see Table 1). For each HAGB, the elastic stiffness tensor $C$ is first rotated in the GB reference frame taking into account the specific GB plane and then in each crystal considering the GB tilt misorientation angle. For the calculations in isotropic elasticity, a Voigt-Reuss-Hill average [108,109] of the cubic elastic constants of the interatomic potentials is performed in order to consider relevant values of $\mu$ (shear modulus) and $\nu$ (Poisson ratio) for further comparisons (see Table 3).

Besides, the FFT spatial grid is regular with a constant step size of approximately $0.2$ Å, while the atoms in the atomistic models are not regularly spaced with an average spacing of $3$ Å. Therefore, in order to be able to compare the two approaches at the exact same positions, an
<figure><img src="image_12.png" /><figcaption>**Fig. 9.** Comparison of the contour plots of the trace of the stress tensor $\sigma_{kk}$ and shear stress $\sigma_{xy}$ for a $\Sigma$149 tilt grain boundary modeled using DSUM, obtained via FDDM-FFT for isotropic and anisotropic Ag. Results were calculated with 1024 $\times$ 1024 pixels.</figcaption></figure>

<table><caption>**Table 3**

Shear modulus $\mu$ and Poisson ratio $\nu$ obtained from the Voigt-Reuss-Hill average [108,109] of the elastic constants associated to the EAM potentials used in the atomistic simulations.</caption><thead><tr><td>System</td><td>μ (GPa)</td><td>ν</td></tr></thead><tbody><tr><td>Al</td><td>30.9</td><td>0.33</td></tr><tr><td>Ag</td><td>29.5</td><td>0.37</td></tr></tbody></table>

interpolation of both atomistic and FFT quantities onto a regular grid should be performed. This new grid is built from the atoms surrounding a given node for atomistic results, and from the different point values surrounding a given node for FFT results. To achieve this, a Gaussian kernel is used as follows [43,110]:

$$ \langle S \rangle = \frac{\Sigma_i (S_i \times W_i)}{\Sigma_i W_i}, \\ W_i = \frac{1}{2\pi\bar{\sigma}^2} \exp\left(-\frac{d_i^2}{2\bar{\sigma}^2}\right), \\ d_i = \sqrt{(x_i - x)^2 + (y_i - y)^2}. \tag{25} $$

In Eq. (25), $\langle S \rangle$ is the weighted averaged (interpolated) quantity, $W_i$ the weighting factors, $\bar{\sigma}$ the standard deviation of the Gaussian function (here the notation $\bar{\sigma}$ is used to make a difference with Eq. (23)), $x_i$ and

$y_i$ the coordinates of the atom i and $x$ and $y$ the coordinates of the considered node in the regular grid. An interpolation cutoff is defined in order to account only for values close to the considered node [43]. This cutoff is set to $3\bar{\sigma}$ while $\bar{\sigma}$ is fixed at $1.1$ Å for Al atomistic systems, $1.2$ Å for Ag atomistic systems and $0.15$ Å for all FFT results. The chosen new spatial step of the common interpolated grid is $0.1$ Å (see Fig. 3).

The computational cost of the FDDM-FFT method is very low. It only represents a few minutes on a single processor for anisotropic elastic calculations with a grid of 1024 $\times$ 1024 pixels, depending on the tolerance used to satisfy local stress equilibrium (see Eq. (A.1) in Appendix A). Here, a good accuracy was reached with a small tolerance of $10^{-6}$ to well approximate the stress fields. In contrast, atomistic simulations require a few tens of hours on a single processor depending on the periodicity of the unit cell to construct the GB. This is clearly the main advantage of the present continuum-based FDDM-FFT method.

## 4. Results: Stress fields of Al and Ag HAGBs

### 4.1. FDDM-FFT calculations vs. analytical solutions in isotropic elasticity

Figs. 4 and 5 present contour plots comparing the isotropic analytical solutions (see Appendix B for details) and the FDDM-FFT results for the trace of the stress tensor, $\sigma_{kk} = \sigma_{xx} + \sigma_{yy} + \sigma_{zz}$, and the shear stress component $\sigma_{xy}$ in the case of Al. These stresses are normalized by K (isotropic elastic bulk modulus) and $2\mu$, respectively, to facilitate direct
<figure><img src="image_13.png" /><figcaption>**Fig. 10.** Contour plot of the interpolated $\sigma_{kk}$ and $\sigma_{xy}$ fields for the Al $\Sigma29$ HAGB, computed from atomistics and FDDM-FFT calculations with the DSUM model.</figcaption></figure>

comparison of the stress fields. For isotropic homogeneous elasticity, the stress polarization tensor $\tau$ reduces to $\tau = C^0 : \varepsilon^{e,\perp}$ (with $C^0$ given by $K$ and $\mu$), which does not need any iteration scheme [34].

Fig. 4 reports the stress contour plots for $\Sigma29$ with structural units $|BC.BC|$ and a misorientation angle $\theta = 46.40^{\circ}$ while Fig. 5 reports the stress contour plots for $\Sigma149$ with structural units $|AABABAB.AABA BAB|$ and a misorientation angle $\theta = 20.02^{\circ}$. The values of the standard deviations, $\sigma$, used in the Gaussian distributions of Eq. (23) were fitted in order to have FFT results as close as possible to the analytical solutions. As can be observed in Figs. 4 and 5, strong agreement between both methods can be achieved without difficulty. Let us note that the calibrated $\sigma$ values are different from one GB to the other with $\sigma = 0.1d$ for the $\Sigma29$ HAGB and $\sigma = 0.6d$ for the $\Sigma149$ HAGB.

The results shown in Figs. 4 and 5 demonstrate that the proposed FDDM-FFT method, which incorporates a Gaussian disclination density distribution, provides highly accurate predictions for stress components when compared to the analytical solutions. It is particularly highlighted by the consistency observed in stress components across the contour plots. It also underscores the efficiency of the FDDM-FFT method for capturing stress distributions of HAGBs. This close match supports its use for detailed stress analysis in isotropic materials. It further demonstrates the robustness of the numerical framework, establishing it as a powerful tool for investigating stress fields with varying disclination characteristics in terms of spatial distribution and strength.

## 4.2. FDDM-FFT predictions of stress fields of HAGBs for Al and Ag: anisotropic vs. isotropic elasticity calculations

In this subsection, we compare the results provided by the FDDM-FFT method considering either isotropic or anisotropic elasticity for two different FFC materials, specifically aluminum (Al) and silver (Ag). The isotropic elastic properties of the materials are given in Table 3 while the anisotropic ones are provided in Table 1. On both sides of the GB, the elastic stiffness tensors of left crystal (denoted $C^{(1)}$) and right crystal (denoted $C^{(2)}$) are rotated by $-\theta/2$ and $\theta/2$ respectively, where $\theta$ is the GB tilt misorientation angle. Therefore, the elastic stiffness tensor is heterogeneous and requires an iterative approach to be determined within the FFT-based method (see Appendix A), with an initial homogeneous medium stiffness $C^0 = 0.5(C^{(1)} + C^{(2)})$ [4]. The analysis focuses on the same stress fields as before ($\sigma_{kk}$ and $\sigma_{xy}$), for both representative HAGBs, namely $\Sigma29$ and $\Sigma149$. Here, since anisotropic elasticity is considered, the stress fields are not normalized by isotropic constants anymore. Contour plots are provided to visualize the stress distributions across these materials, highlighting the differences in stress profiles between the isotropic and anisotropic cases. The FDDM-FFT method captures these stress components under varying elastic properties. Same parameters for Gaussian distributions as for isotropic elasticity were used in the case of FFT-based anisotropic calculations.

Figs. 6 and 7 illustrate the comparison between isotropic and anisotropic FDDM-FFT results for Al and Ag $\Sigma29$. For this HAGB, stress
<figure><img src="image_14.png" /><figcaption>Fig. 11. Contour plot of the interpolated $\sigma_{kk}$ and $\sigma_{xy}$ fields for the Al $\Sigma149$ HAGB, computed from atomistics and FDDM-FFT calculations with the DSUM model.</figcaption></figure>

equilibrium was reached after 5 iterations for Al and after 14 iterations for Ag (using Eq. (A.1) of Appendix B as convergence criterion).

Figs. 8 and 9 report the same types of comparison with the FDDM-FFT method for $\Sigma149$. For this HAGB, stress equilibrium was reached after 5 iterations for Al and after 10 iterations for Ag (using Eq. (A.1) of Appendix B as convergence criterion).

Aluminum (Al) can be considered nearly isotropic with a Zener ratio of 1.22. Consistently, Figs. 6 and 8 show minimal differences between the isotropic and anisotropic cases for both $\Sigma29$ and $\Sigma149$ HAGBs. However, silver (Ag) has a much higher Zener ratio of 3.21. Accordingly, the influence of elastic anisotropy becomes more pronounced for both $\Sigma29$ and $\Sigma149$, leading to clear differences between isotropic and anisotropic behaviors, as can be analyzed from Figs. 7 and 9. While contour plots for anisotropic Al are still rounded as for the isotropic case, the contour plots for Ag are much sharper. In the next subsection, the FDDM-FFT calculations with heterogeneous anisotropic elasticity are presented with some comparisons to atomistic results.

### 4.3. FDDM-FFT calculations vs. atomistic simulations

The anisotropic FDDM-FFT stress fields of $\Sigma29(520)[001]46.40^{\circ}$ and $\Sigma149(1070)[001]20.02^{\circ}$, obtained with the DSUM and presented in Section 4.2, are now compared to virial stress fields evaluated using molecular statics (MS). In both cases, stress fields are interpolated onto a regular grid with a spatial resolution of 0.1 Å, following the methodology described in Section 3.3. Comparisons are performed for both Al and Ag systems (see Table 1). Figs. 10 and 11 first show interpolated 2D contour plots of the trace of the stress tensor $\sigma_{kk}$ and

the shear stress $\sigma_{xy}$, obtained via MS computation and DSUM for Al. The contours of $\sigma_{xy}$ derived from DSUM closely align those obtained by atomistic simulations, for both $\Sigma29$ and $\Sigma149$ HAGBs. However, discrepancies arise from the Gaussian representation of core defects in DSUM. They are more pronounced for $\sigma_{kk}$, and particularly in the $\Sigma29$ GB, where atomistic fields exhibits a different core shape extension. In contrast, the contour plots for $\Sigma149$ GB exhibit stronger agreement between both models.

Despite these differences, the interpolated atomistic and DSUM stress field variations exhibit the same periodicity for all grain boundaries, as illustrated from stress profiles along y-axis at 3 Å from GB in Fig. 12. Especially for the $\Sigma149$ GB, FDDM-FFT calculations with DSUM model captures very well the variations in stress profiles for both Al and Ag. However, differences in magnitude do persist. This is due to the fact that the Gaussian parameter $\sigma$ was only calibrated using the analytical solution in isotropic elasticity (see Section 4.1).

To improve the agreement between both methods, especially regarding the stress field magnitudes, the parameter $\sigma$ of the Gaussian distribution function could be calibrated again, but now based on these atomistic results. Figs. 13 and 14 show 2D contour plots of the stress fields for Ag GBs. The Zener ratio being significantly higher in Ag than in Al (see Table 1), it is reflected in the stress field distribution (see Figs. 10, 13, 11 and 14). Compared to Al, the stress fields in Ag are more extended and exhibit a less circular shape around the defect core. This behavior, indicating a higher elastic anisotropy effect, is well captured by the DSUM. By comparing Figs. 6 and 8, Figs. 7 and 9, and Fig. 13 (or Fig. 14), it is highlighted that the effect of elastic anisotropy leads to a pronounced breaking of symmetry in the plane normal to the y-axis, which is in agreement with the MS results. Fig. 15, which displays
<figure><img src="image_15.png" /><figcaption>**Fig. 12.** Comparison of the interpolated profiles of $\sigma_{kk}$ and $\sigma_{xy}$ along y-axis at $3\AA$ from GB (see inset) for the $\Sigma29$ and $\Sigma149$ HAGBs in Al computed from atomistics and FDDM-FFT calculations with the DSUM model.</figcaption></figure>

stress profiles along y-axis at $3\AA$ from GB, further demonstrates that the DSUM effectively captures the anisotropic elastic behavior of Ag. As previously, calibrating again the Gaussian parameter $\sigma$ using atomistic data would be necessary to achieve a better match with respect to stress magnitude profiles. This parametric study is left for future. Figs. 12 and 15 display a better agreement for the shear stress $\sigma_{xy}$ compared to the stress tensor trace $\sigma_{kk}$ (which contains the three normal components: $\sigma_{xx}$, $\sigma_{yy}$ and $\sigma_{zz}$). This difference is due to the cumulative errors of the sum of these three components for $\sigma_{kk}$, while $\sigma_{xy}$ contains only one component leading to a smaller error.

## 5. Conclusions

In conclusion, this work reveals the effect of anisotropic elasticity on the elastic fields of two different High Angle Grain Boundaries (HAGBs) considering two different FCC systems, with low (Al) and medium (Ag) Zener ratios. Continuum-based theoretical predictions were performed using a FFT-based method in the framework of dislocation and disclination mechanics (FDDM-FFT), extending the work of Berbenni and Taupin [37] to HAGBs, with a description of HAGBs in terms of disclination structural units (DSUM) constructed from atomistics.

This was done by calibrating the strengths of disclination dipoles in the DSUM, as well as their spread with Gaussian distributions, in comparisons with analytical isotropic elastic solutions.

Then, the main challenge was to compare the elastic fields obtained from the FDDM-FFT method in heterogeneous anisotropic elasticity, with atomistic simulation results onto the same spatial grid. This was achieved using an interpolation method able to extract continuum stress fields from atomic virial stresses. As a result, it was observed that the present FDDM-FFT method based on the DSUM model captures well the main features of the anisotropic elastic fields obtained by atomistics for such HAGBs (shape of contour maps, stress profiles, periodicity). This work represents an extension of the anisotropic elastic calculations of Petrazoller et al. [43], who used the Stroh formalism with dislocation models for Low Angle Grain Boundaries (LAGBs).

As an improvement of this method, a more refined description of the defect core extension will be studied by calibrating the spatial distribution of disclination densities directly from interpolated fields extracted from atomistic data.

As a perspective to this work, the present FDDM-FFT anisotropic framework will be applied to model solute segregation at HAGBs in order to extend our previous contributions in this research area [38,
<figure><img src="image_16.png" /><figcaption>Fig. 13. Contour plot of the interpolated $\sigma_{kk}$ and $\sigma_{xy}$ fields for the Ag Σ29 HAGB, computed from atomistics and FDDM-FFT calculations with the DSUM model.</figcaption></figure>

[43,69]. Another future work will focus on analyzing stress driven GB motion like shear-coupled GB migration [67,111] in the framework of the present FDDM-FFT method and studying the effect of elastic anisotropy on this mechanism for HAGBs.

### CRediT authorship contribution statement

**Abdallah Wazne:** Writing – original draft, Visualization, Software, Investigation, Formal analysis, Data curation, Conceptualization. **Joé Petrazoller:** Writing – review & editing, Visualization, Software, Investigation, Formal analysis, Data curation, Conceptualization. **Julien Guénolé:** Writing – review & editing, Validation, Supervision, Methodology, Funding acquisition, Conceptualization. **Thiebaud Richeton:** Writing – review & editing, Validation, Supervision, Project administration, Methodology, Funding acquisition, Conceptualization. **Stéphane Berbenni:** Writing – review & editing, Writing – original draft, Validation, Supervision, Software, Methodology, Conceptualization.

### Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

The authors are grateful to the French National Research Agency (ANR) for financial support under contracts ANR-22-CE92-0058-01

(SILA) and ANR-21-CE08-0001 (ATOUUM). The authors also acknowledge the use of the Cassiopee Arts et Métiers Institute of Technology HPC Center made available for conducting the research reported in this paper.

## Appendix A. FFT-based algorithm for elasto-static field dislocation and disclination mechanics (FDDM)

The direct and inverse Fourier transforms are performed using the Fast Fourier Transform (FFT) algorithm [112]. The unit cell $V$ is assumed to have spatial dimensions $T_1, T_2$, and $T_3$ along the $e_1, e_2$, and $e_3$ axes, respectively. It is discretized on a regular rectangular grid with $N_1 \times N_2 \times N_3$ voxels, with position vector $\mathbf{x} = (i_1\delta_1, i_2\delta_2, i_3\delta_3)$, where $i_1 = 0 \rightarrow N_1-1$, $i_2 = 0 \rightarrow N_2-1$, $i_3 = 0 \rightarrow N_3-1$, and $\delta_1, \delta_2, \delta_3$ represent the voxel sizes in the $e_1, e_2$, and $e_3$ directions, respectively. In the case of three-dimensional (3D) problems, the total number of FFT grid points is $N_{tot} = N_1 \times N_2 \times N_3$, with $d = \delta_1 = \delta_2 = \delta_3$. The discrete Fourier transform (DFT) of a given spatial function $f$ is denoted by $\hat{f} = \text{FFT}(f)$, and its inverse is given by $f = \text{FFT}^{-1}(\hat{f})$. For the applications presented in Section 4, we will focus on two-dimensional (2D) cases involving bulk and interfacial line defects, where $N_{tot} = N_1 \times N_2$ and $d = \delta_1 = \delta_2$, to illustrate the method.

The algorithm proceeds as follows: After specifying the disclination density $\theta(\mathbf{x})$ and the dislocation density $\alpha(\mathbf{x})$ in real space, the initialization begins by calculating $\hat{\theta}(\xi)$ and $\hat{\alpha}(\xi)$ using a direct FFT (step 1). Next, Eqs. (20) are applied to compute $\hat{\kappa}^{e,\perp}(\xi)$ and $\hat{\kappa}^{e,\perp}(\xi)$ in Fourier space (step 2), followed by the computation of $\epsilon^{e,\perp}(\mathbf{x})$ in real space using the inverse FFT (step 3). The initial compatible elastic strain $\epsilon_0^{e,\parallel}$
<figure><img src="image_17.png" /><figcaption>Fig. 14. Contour plot of the interpolated $\sigma_{kk}$ and $\sigma_{xy}$ fields for the Ag $\Sigma$149 HAGB, computed from atomistics and FDDM-FFT calculations with the DSUM model.</figcaption></figure>

**Algorithm 1** FFT-based algorithm for elasto-static field equations of FDDM

**Initialization:** $\theta(\mathbf{x})$ and $\alpha(\mathbf{x})$ are prescribed

1. $\hat{\theta}(\xi) \leftarrow \text{FFT}(\theta)$ and $\hat{\alpha}(\xi) \leftarrow \text{FFT}(\alpha)$

2. Solve Poisson equations to obtain $\hat{\varepsilon}^{e,\perp}(\xi)$ and $\hat{\kappa}^{e,\perp}(\xi)$:

$$ (i) \quad \hat{\kappa}_{ji}^{e,\perp}(\xi) = i \frac{\xi_k}{\xi^2} e_{jkl} \hat{\theta}_{ll}(\xi) \quad \forall \xi \neq 0, \quad \hat{\kappa}_{ji}^{e,\perp}(0) = 0 $$

$$ (ii) \quad \hat{\varepsilon}_{(lr)}^{e,\perp}(\xi) = i \frac{\xi^q}{\xi^2} e_{rql} \left( \hat{\alpha}_{il}(\xi) - \delta_{il} \hat{\kappa}_{pp}^{e,\perp}(\xi) + \hat{\kappa}_{il}^t e^{e,\perp}(\xi) \right)_{(lr)} \quad \forall \xi \neq 0 $$

3. $\varepsilon^{e,\perp}(\mathbf{x}) \leftarrow \text{FFT}^{-1}(\hat{\varepsilon}^{e,\perp})$

4. $\varepsilon_0^{e,\parallel} \leftarrow \langle \varepsilon^{e,\parallel} \rangle_0 = C^{0-1} : \bar{\sigma}$

5. $\sigma_0(\mathbf{x}) \leftarrow C(\mathbf{x}) : (\varepsilon_0^{e,\parallel} + \varepsilon^{e,\perp}(\mathbf{x}))$

**Iterate:** $n+1$ ($\varepsilon_n^{e,\parallel}(\mathbf{x})$ and $\sigma_n(\mathbf{x})$ being known)

6. $\hat{\sigma}_n(\xi) = \text{FFT}(\sigma_n)$

7. Convergence test based on Eq. (A.1)

8. $\hat{\tau}_n(\xi) \leftarrow \hat{\sigma}_n(\xi) - C^0 : \varepsilon_n^{e,\parallel}(\xi)$

9. $\hat{\varepsilon}_{n+1}^{e,\parallel}(\xi) \leftarrow -\hat{\Gamma}^0(\xi) : \hat{\tau}_n(\xi) \quad \forall \xi \neq 0$ and $\hat{\varepsilon}_{n+1}^{e,\parallel}(0) \leftarrow \langle \varepsilon^{e,\parallel} \rangle_n$

10. $\varepsilon_{n+1}^{e,\parallel}(\mathbf{x}) \leftarrow \text{FFT}^{-1}(\hat{\varepsilon}_{n+1}^{e,\parallel}(\xi))$

11. $\sigma_{n+1}(\mathbf{x}) \leftarrow C(\mathbf{x}) : (\varepsilon_{n+1}^{e,\parallel}(\mathbf{x}) + \varepsilon^{e,\perp}(\mathbf{x}))$

In the global iterative loop, at iteration $(n+1)$, the stress field from the previous iteration, denoted $\sigma_n(\mathbf{x})$, is transformed into Fourier space using a direct FFT to obtain $\hat{\sigma}_n(\xi)$ (step 6). This transformed stress field is then used to check the convergence criterion based on stress equilibrium in Fourier space (step 7). If convergence is achieved, the iteration is terminated; otherwise, the polarization stress tensor $\hat{\tau}_n(\xi)$ is computed in Fourier space by calculating the FFT of $C^0 : \varepsilon_n^{e,\parallel}$ and subtracting it from $\hat{\sigma}_n(\xi)$ (step 8). The compatible elastic strain is then computed in Fourier space at step 9 using Eq. (21), resulting in $\hat{\varepsilon}_{n+1}^{e,\parallel}(\xi)$ ($\forall \xi \neq 0$). For $\xi = 0$, $\hat{\varepsilon}_{n+1}^{e,\parallel}(0)$ is set to the average value $\langle \varepsilon^{e,\parallel} \rangle$ from iteration $(n)$, considering the traction boundary conditions as in [35]. The compatible elastic strain is then transformed back to real space via inverse FFT (step 10). This updated strain is used to update the stress field in real space (step 11). The iterative loop continues with the updated strain and stress fields until convergence is achieved, as described in Eq. (A.1). Following [4], the convergence criterion for stress equilibrium expressed in Fourier space is applied at step 7 and is expressed as:

$$ e_n = \frac{\left( \frac{1}{N} \sum_{i=1}^{N} |\xi_i \cdot \hat{\sigma}_n(\xi_i)|^2 \right)^{1/2}}{\| \hat{\sigma}_n(0) \|_F} \le \epsilon, \quad (A.1) $$

where $\|\cdot\|_F$ represents the Frobenius norm for second order tensors, $|\cdot|$ is the absolute value, and $\epsilon_n$ is the error at iteration $n$. Convergence is achieved when the error falls below a specified precision $\epsilon$ (set to $10^{-6}$ in our simulations). Once convergence is reached, the elastic strain (including both incompatible and compatible components) and the stress field are determined in real space.

is taken as the homogeneous elastic solution $C^{0-1} : \bar{\sigma}$ (step 4). The initialization concludes with the calculation of the initial stress field $\sigma_0(\mathbf{x})$, which is given by $C(\mathbf{x}) : (\varepsilon_0^{e,\parallel} + \varepsilon^{e,\perp}(\mathbf{x}))$ (step 5).
<figure><img src="image_18.png" /><figcaption>Fig. 15. Comparison of the interpolated profiles of $\sigma_{kk}$ and $\sigma_{xy}$ along y-axis at 3Å from GB (see inset) for the Σ29 and Σ149 HAGBs in Ag computed from atomistics and FDDM-FFT calculations with the DSUM model.</figcaption></figure>

An important issue in FFT-based methods is to predict accurate local stress fields near material discontinuities, such as defects and inhomogeneities, avoiding numerical oscillations that arise from the classic FFT approximation. Therefore, finite difference (FD) schemes for computing first- and second-order spatial derivatives were introduced alongside the discrete Fourier transform (DFT). For the present simulations, a centered FD scheme for second derivatives together with DFT as in Berbenni et al. [34], was employed to solve the Poisson equations in the Fourier space (at step 2), while a rotated centered FD scheme as proposed by Willot [113] was used to compute the Green operator in the Fourier space for the Lippmann-Schwinger equation (at step 9). The combination of both schemes were seen to give accurate results in the case of field dislocation mechanics calculations for strongly heterogeneous elastic materials, as for instance a dislocation in a matrix interacting with a void within a matrix phase [35].

## Appendix B. Analytical formulae for stress field of disclination dipole wall in isotropic elasticity

Building on the works by Hurtado et al. [75], and Dingreville and Berbenni [69], the stress field, $\sigma_{ij}^*(x, y)$, of an infinite disclination wall

in the yz-plane, made up of disclination dipoles lying parallel to the z-axis with same strength $\omega$, same arm length $0.5L_n$ and spaced at a periodic distance $H$, can be expressed in linear isotropic elasticity as:

$$ \frac{1}{2} (\sigma_{xx}^{\ddagger} + \sigma_{yy}^{\ddagger}) = \frac{1}{2} D_0 \omega \ln \frac{\cosh \alpha - \cos (\beta - \lambda_n)}{\cosh \alpha - \cos (\beta + \lambda_n)} \quad (B.1) $$

$$ \frac{1}{2} (\sigma_{xx}^{\ddagger} - \sigma_{yy}^{\ddagger}) = \frac{1}{2} D_0 \omega \alpha \left[ \frac{\sinh \alpha}{\cosh \alpha - \cos (\beta + \lambda_n)} - \frac{\sinh \alpha}{\cosh \alpha - \cos (\beta - \lambda_n)} \right] \quad (B.2) $$

$$ \sigma_{xy}^{\ddagger} = \frac{1}{2} D_0 \omega \alpha \left[ \frac{\sin (\beta + \lambda_n)}{\cosh \alpha - \cos (\beta + \lambda_n)} - \frac{\sin (\beta - \lambda_n)}{\cosh \alpha - \cos (\beta - \lambda_n)} \right] \quad (B.3) $$

with $\alpha = 2\pi x/H$, $\beta = 2\pi y/H$ and $\lambda_n = \pi L_n/H$. $D_0 = \mu/2\pi(1-\nu)$ where $\mu$ is the elastic shear modulus and $\nu$ the elastic Poisson ratio.

## Data availability

Data will be made available on request.
References

[1] M. Schneider, A review of nonlinear FFT-based computational homogenization methods, Acta Mech. 232 (6) (2021) 2051-2100, http://dx.doi.org/10.1007/s00707-021-02962-1.

[2] S. Lucarini, M.V. Upadhyay, J. Segurado, FFT based approaches in micromechanics: fundamentals, methods and applications, Modelling Simul. Mater. Sci. Eng. 30 (2) (2021) 023002, http://dx.doi.org/10.1088/1361-651X/ac34e1.

[3] H. Moulinec, P. Suquet, A fast numerical method for computing the linear and non linear properties of composites, C. R. Acad. Sci. Paris II 318 (1994) 1417-1423.

[4] H. Moulinec, P. Suquet, A numerical method for computing the overall response of nonlinear composites with complex microstructure, Comput. Methods Appl. Mech. Engrg. 157 (1998) 69-94, http://dx.doi.org/10.1016/S0045-7825(97) 00218-1.

[5] R. Lebensohn, N-site modeling of a 3D viscoplastic polycrystal using fast Fourier transform, Acta Mater. 49 (2001) 2723-2737, http://dx.doi.org/10.1016/S1359-6454(01)00172-0.

[6] R. Lebensohn, R. Brenner, O. Castelnau, A. Rollett, Orientation image-based micromechanical modelling of subgrain texture evolution in polycrystalline copper, Acta Mater. 56 (2008) 3914-3926, http://dx.doi.org/10.1016/j.actamat. 2008.04.016.

[7] R.A. Lebensohn, A.K. Kanjarla, P. Eisenlohr, An elasto-viscoplastic formulation based on fast Fourier transforms for the prediction of micromechanical fields in polycrystalline materials, Int. J. Plast. 32-33 (2012) 59-69, http://dx.doi.org/10.1016/j.ijplas.2011.12.005.

[8] P. Eisenlohr, M. Diehl, R.A. Lebensohn, F. Roters, A spectral method solution to crystal elasto-viscoplasticity at finite strains, Int. J. Plast. 46 (2013) 37-53, http://dx.doi.org/10.1016/j.ijplas.2012.09.012.

[9] R.A. Lebensohn, A.D. Rollett, Spectral methods for full-field micromechanical modelling of polycrystalline materials, Comput. Mater. Sci. 173 (2020) 109336, http://dx.doi.org/10.1016/j.commatsci.2019.109336.

[10] A. Prakash, R. Lebensohn, Simulation of micromechanical behavior of polycrystals: finite elements versus fast Fourier transforms, Model. Simul. Mater. Sci. Eng. 17 (2009) 64010-64016, http://dx.doi.org/10.1088/0965-0393/17/ 6/064010.

[11] B. Liu, D. Raabe, F. Roters, P. Eisenlohr, R. Lebensohn, Comparison of finite element and fast Fourier transform crystal plasticity solvers for texture prediction, Modelling Simul. Mater. Sci. Eng. 18 (2009) 085005, http://dx.doi.org/10.1088/0965-0393/18/8/085005.

[12] B.S. Anglin, R.A. Lebensohn, A.D. Rollett, Validation of a numerical method based on fast Fourier transforms for heterogeneous thermoelastic materials by comparison with analytical solutions, Comput. Mater. Sci. 87 (2014) 209-217, http://dx.doi.org/10.1016/j.commatsci.2014.02.027.

[13] S.P. Donegan, A.D. Rollett, Simulation of residual stress and elastic energy density in thermal barrier coatings using fast Fourier transforms, Acta Mater. 96 (2015) 212-228, http://dx.doi.org/10.1016/j.actamat.2015.06.019.

[14] K.S. Eloh, A. Jacques, S. Berbenni, Development of a new consistent discrete green operator for FFT-based methods to solve heterogeneous problems with eigenstrains, Int. J. Plast. 116 (2019) 1-23, http://dx.doi.org/10.1016/j.ijplas. 2018.10.011.

[15] R. Santos-Güemes, C.J. Ortiz, J. Segurado, An FFT based approach to account for elastic interactions in OkMC: Application to dislocation loops in iron, J. Nucl. Mater. 594 (2024) 155020, http://dx.doi.org/10.1016/j.jnucmat.2024. 155020.

[16] R. Lebensohn, J. Escobedo, E. Cerreta, D. Dennis-Koller, C. Bronkhorst, J. Bingert, Modeling void growth in polycrystalline materials, Acta Mater. 61 (2013) 6918-6932, http://dx.doi.org/10.1016/j.actamat.2013.08.004.

[17] R.A. Lebensohn, A. Needleman, Numerical implementation of non-local polycrystal plasticity using fast Fourier transforms, J. Mech. Phys. Solids 97 (2016) 333-351, http://dx.doi.org/10.1016/j.jmps.2016.03.023.

[18] P.G. Christodoulou, R.A. Lebensohn, I.J. Beyerlein, Study of the interplay between lower-order and higher-order energetic strain-gradient effects in polycrystal plasticity, J. Mech. Phys. Solids 164 (2022) 104906, http://dx.doi.org/ 10.1016/j.jmps.2022.104906.

[19] M. Zecevic, R.A. Lebensohn, L. Capolungo, Non-local large-strain FFT-based formulation and its application to interface-dominated plasticity of nano-metallic laminates, J. Mech. Phys. Solids 173 (2023) 105187, http://dx.doi.org/10.1016/j.jmps.2022.105187.

[20] M. Zecevic, R.A. Lebensohn, L. Capolungo, Generalized grain boundary constitutive description implemented in a strain-gradient large-strain FFT-based formulation: Application to nano-metallic laminates, J. Mech. Phys. Solids 193 (2024) 105859, http://dx.doi.org/10.1016/j.jmps.2024.105859.

[21] N.M. Francis, F. Pourahmadian, R.A. Lebensohn, R. Dingreville, A fast Fourier transform-based solver for elastic micropolar composites, Comput. Methods Appl. Mech. Engrg. 418 (2024) 116510, http://dx.doi.org/10.1016/j.cma.2023. 116510.

[22] N.M. Francis, R.A. Lebensohn, F. Pourahmadian, R. Dingreville, Micropolar elastoplasticity using a fast Fourier transform-based solver, Internat. J. Numer. Methods Engrg. 126 (1) (2025) e7651, http://dx.doi.org/10.1002/nme.7651.

[23] Y. Chen, D. Vasiukov, L. Gélébart, C.H. Park, A FFT solver for variational phase-field modeling of brittle fracture, Comput. Methods Appl. Mech. Engrg. 349 (1967) 169-190, http://dx.doi.org/10.1016/j.cma.2019.02.017.

[24] F. Ernesti, M. Schneider, T. Böhlke, Fast implicit solvers for phase-field fracture problems on heterogeneous microstructures, Comput. Methods Appl. Mech. Engrg. 363 (2020) 112793, http://dx.doi.org/10.1016/j.cma.2019.112793.

[25] S. Lucarini, J. Segurado, On the accuracy of spectral solvers for micromechanics based fatigue modeling, Comput. Mech. 63 (2019) 365-382, http://dx.doi.org/ 10.1007/s00466-018-1598-1.

[26] K.S. Djaka, S. Berbenni, V. Taupin, R.A. Lebensohn, A FFT-based numerical implementation of mesoscale field dislocation mechanics: Application to two-phase laminates, Int. J. Solids Struct. 184 (2020) 136-152, http://dx.doi.org/ 10.1016/j.ijsolstr.2018.12.027.

[27] S. Berbenni, V. Taupin, R.A. Lebensohn, A fast Fourier transform-based mesoscale field dislocation mechanics study of grain size effects and reversible plasticity in polycrystals, J. Mech. Phys. Solids 135 (2020) 103808, http://dx.doi.org/10.1016/j.jmps.2019.103808.

[28] J. Genée, N. Gey, F. Bonnet, R.A. Lebensohn, S. Berbenni, Experimental and numerical investigation of key microstructural features influencing the localization of plastic deformation in Fe-TiB2 metal matrix composite, J. Mater. Sci. 56 (19) (2021) 11278-11297, http://dx.doi.org/10.1007/s10853-021-06017-7.

[29] S. Berbenni, R.A. Lebensohn, A numerical study of reversible plasticity using continuum dislocation mechanics, C. R. Phys. 22 (S3) (2021) 295-312, http://dx.doi.org/10.5802/crphys.v54.

[30] S. Berbenni, V. Taupin, R.A. Lebensohn, A FFT-based mesoscale continuum dislocation mechanics with defect energy: Applications to composites and polycrystals, Eur. J. Mech. A Solids 111 (2025) 105548, http://dx.doi.org/10. 1016/j.euromechsol.2024.105548.

[31] J.T. Graham, A.D. Rollett, R. LeSar, Fast Fourier transform discrete dislocation dynamics, Model. Simul. Mater. Sci. Eng. 8 (2016) 085005, http://dx.doi.org/ 10.1088/0965-0393/24/8/085005.

[32] N. Bertin, L. Capolungo, A FFT-based formulation for discrete dislocation dynamics in heterogeneous media, J. Comput. Phys. 355 (2018) 366-384, http://dx.doi.org/10.1016/j.jcp.2017.11.020.

[33] R. Brenner, A.J. Beaudoin, P. Suquet, A. Acharya, Numerical implementation of static field dislocation mechanics theory for periodic media, Phil. Mag. (2014) http://dx.doi.org/10.1080/14786435.2014.896081.

[34] S. Berbenni, V. Taupin, K.S. Djaka, C. Fressengeas, A numerical spectral approach for solving elasto-static field dislocation and G-disclination mechanics, Int. J. Solids Struct. 51 (2014) 4157-4175, http://dx.doi.org/10.1016/j.ijsolstr. 2014.08.009.

[35] K.S. Djaka, A. Villani, V. Taupin, L. Capolungo, S. Berbenni, Field dislocation mechanics for heterogeneous elastic materials: A numerical spectral approach, Comput. Methods Appl. Mech. Engrg. 315 (2017) 921-942, http://dx.doi.org/ 10.1016/j.cma.2016.11.036.

[36] S. Berbenni, V. Taupin, C. Fressengeas, L. Capolungo, A fast Fourier transform-based approach for generalized disclination mechanics within a couple stress theory, in: H. Altenbach, S. Forest (Eds.), Generalized Continua as Models for Classical and Advanced Materials, in: Advanced Structured Materials, Springer International Publishing Switzerland, 2016, pp. 47-75.

[37] S. Berbenni, V. Taupin, Fast Fourier transform-based micromechanics of interfacial line defects in crystalline materials, J. Micromech. Mol. Phys. 1840007 (2018) http://dx.doi.org/10.1142/S2424913018400076.

[38] C. Hu, S. Berbenni, D.L. Medlin, R. Dingreville, Discontinuous segregation patterning across disconnections, Acta Mater. 246 (2023) 118724, http://dx. doi.org/10.1016/j.actamat.2023.118724.

[39] A.P. Sutton, R.W. Balluffi, Interfaces in Crystalline Solids, Oxford University Press, Oxford (UK), 1995.

[40] P.M. Anderson, J.P. Hirth, J. Lothe, Theory of Dislocations, third ed., Cambridge University Press, 2017.

[41] W.T. Read, W. Shockley, Dislocation models of crystal grain boundaries, Phys. Rev. 78 (1950) 275-289, http://dx.doi.org/10.1103/PhysRev.78.275.

[42] W.T. Read, The motion of a grain boundary, in: Dislocations in Crystals, Mc Graw-Hill, New York, 1953, (Chapter 14).

[43] J. Petrazoller, J. Guénolé, S. Berbenni, T. Richeton, On the effect of elastic anisotropy and polarizability on solute segregation at low-angle grain boundaries, Comput. Mater. Sci. 249 (2025) 113642, http://dx.doi.org/10.1016/j. commatsci.2024.113642.

[44] A.N. Stroh, Dislocations and cracks in anisotropic elasticity, Phil. Mag. 3 (1958) 625-646, http://dx.doi.org/10.1080/14786435808565804.

[45] D.M. Heyes, Pressure tensor of partial-charge and point-dipole lattices with bulk and surface geometries, Phys. Rev. B 49 (2) (1994) 755-764, http://dx.doi.org/ 10.1103/PhysRevB.49.755.

[46] A.P. Thompson, S.J. Plimpton, W. Mattson, General formulation of pressure and stress tensor for arbitrary many-body interaction potentials under periodic boundary conditions, J. Chem. Phys. 131 (15) (2009) 154107, http://dx.doi. org/10.1063/1.3245303, Publisher: AIP Publishing.

[47] A. Stukowski, A. Arsenlis, On the elastic-plastic decomposition of crystal deformation at the atomic scale, Modelling Simul. Mater. Sci. Eng. 20 (3) (2012) 035012, http://dx.doi.org/10.1088/0965-0393/20/3/035012, Publisher: IOP Publishing.
[48] P.M. Bollmann, J.P. Hirth, J. Lothe, Theory of Dislocations, third ed., Springer, Berlin, 2017.

[49] F.C. Frank, Symposium on the Plastic Deformation of Crystalline Solids, Carnegie Institute of Technology, Pittsburgh, 1950, pp. 150–154.

[50] B.A. Bilby, Bristol Conference Report on Defects in Crystalline Solids, The Physical Society of London, 1955, pp. 124–133.

[51] R. Bullough, B.A. Bilby, Continuous distributions of dislocations: surface dislocations and the crystallography of martensitic transformations, Proc. Phys. Soc. B 69 (1956) 1276–1286, http://dx.doi.org/10.1088/0370-1301/69/12/311.

[52] J. Hirth, R. Pond, R. Hoagland, X.-Y. Liu, R. Wang, Interface defects, reference spaces and the Frank–Bilby equation, Prog. Mater. Sci. 58 (5) (2013) 749–823, http://dx.doi.org/10.1016/j.pmatsci.2012.10.002.

[53] T. Mura, Micromechanics of Defects in Solids, Kluwer Academic Publishers, Dordrecht, The Netherlands, 1987.

[54] G. Bishop, B. Chalmers, A coincidence — Ledge — Dislocation description of grain boundaries, Scr. Metall. 2 (2) (1968) 133–139, http://dx.doi.org/10.1016/0036-9748(68)90085-9.

[55] A.P. Sutton, V. Vitek, On the structure of tilt grain boundaries in cubic metals I. Symmetrical tilt boundaries, Philos. Trans. R. Soc. Lond. Ser. A Math. Phys. Sci. 309 (1983) 1–36, http://dx.doi.org/10.1098/rsta.1983.0020.

[56] A.P. Sutton, V. Vitek, On the structure of tilt grain boundaries in cubic metals II. Asymmetrical tilt boundaries, Philos. Trans. R. Soc. Lond. Ser. A Math. Phys. Sci. 309 (1983) 37–54, http://dx.doi.org/10.1098/rsta.1983.0021.

[57] A.P. Sutton, V. Vitek, On the structure of tilt grain boundaries in cubic metals III. Generalizations of the structural study and implications for the properties of grain boundaries, Philos. Trans. R. Soc. Lond. Ser. A Math. Phys. Sci. 309 (1983) 55–68, http://dx.doi.org/10.1098/rsta.1983.0022.

[58] J. Han, V. Vitek, D.J. Srolovitz, Grain-boundary metastability and its statistical properties, Acta Mater. 104 (2016) 259–273, http://dx.doi.org/10.1016/j.actamat.2015.11.035.

[59] J. Han, V. Vitek, D.J. Srolovitz, The grain-boundary structural unit model redux, Acta Mater. 133 (2017) 186–199, http://dx.doi.org/10.1016/j.actamat.2017.05.002, URL https://www.sciencedirect.com/science/article/pii/S1359645417303798.

[60] C. Hartley, Y. Mishin, Characterization and visualization of the lattice misfit associated with dislocation cores, Acta Mater. 53 (5) (2005) 1313–1321, http://dx.doi.org/10.1016/j.actamat.2004.11.027.

[61] G.J. Tucker, J.A. Zimmerman, D.L. McDowell, Continuum metrics for deformation and microrotation from atomistic simulations: application to grain boundaries, Internat. J. Engrg. Sci. 49 (2011) 1424–1434, http://dx.doi.org/10.1016/j.ijengsci.2011.03.019.

[62] P. van Beers, V. Kouznetsova, M. Geers, M. Tschopp, D. McDowell, A multiscale model of grain boundary structure and energy: From atomistics to a continuum description, Acta Mater. 82 (2015) 513–529, http://dx.doi.org/10.1016/j.actamat.2014.08.045.

[63] D.H. Warner, F. Sansoz, J.F. Molinari, Atomistic based continuum investigation of plastic deformation in nanocrystalline copper, Int. J. Plast. 22 (2006) 754–774, http://dx.doi.org/10.1016/j.ijplas.2005.04.014.

[64] J.D. Clayton, D.L. McDowell, D.J. Bammann, Modeling dislocations and dislocations with finite micropolar elastoplasticity, Int. J. Plast. 22 (2006) 210–256, http://dx.doi.org/10.1016/j.ijplas.2004.12.001.

[65] J. Wang, A. Misra, J.P. Hirth, Shear response of twin boundaries in face-centered-cubic metals, Phys. Rev. B 83 (2011) 064106, http://dx.doi.org/10.1103/PhysRevB.83.064106.

[66] M. McPhie, S. Berbenni, M. Cherkaoui, Activation energy for nucleation of partial dislocation from grain boundaries, Comput. Mater. Sci. 62 (2012) 169–174, http://dx.doi.org/10.1016/j.commatsci.2012.05.039.

[67] S. Berbenni, B. Paliwal, M. Cherkaoui, A micromechanics-based model for shear-coupled grain boundary migration in bicrystals, Int. J. Plast. 44 (2013) 68–94, http://dx.doi.org/10.1016/j.ijplas.2012.11.011.

[68] A. Prakash, J. Guénolé, J. Wang, J. Müller, E. Spiecker, M.J. Mills, I. Povstugar, P. Choi, D. Raabe, E. Bitzek, Atom probe informed simulations of dislocation-precipitate interactions reveal the importance of local interface curvature, Acta Mater. 92 (2015) 33–45, http://dx.doi.org/10.1016/j.actamat.2015.03.050, arXiv: 1707.04192 ISBN: 1359-6454.

[69] R. Dingreville, S. Berbenni, On the interaction of solutes with grain boundaries, Acta Mater. 104 (2016) 237–249, http://dx.doi.org/10.1016/j.actamat.2015.11.017.

[70] V.Y. Gertsman, A.A. Nazarov, A.E. Romanov, R.Z. Valiev, V.I. Vladimirov, Dislocation-structural unit model of grain boundaries, Philos. Mag. A 59 (5) (1989) 1113–1118, http://dx.doi.org/10.1080/01418618908209841.

[71] A.A. Nazarov, O.A. Shenderova, D.W. Brenner, On the dislocation-structural unit model of grain boundaries, Mater. Sci. Eng. A 281 (1) (2000) 148–155, http://dx.doi.org/10.1016/S0921-5093(99)00727-3.

[72] A.A. Nazarov, A.E. Romanov, On the average misorientation angle of general tilt boundaries, Phil. Mag. Lett. 60 (5) (1989) 187–193, http://dx.doi.org/10.1080/09500838908206456.

[73] J.C.M. Li, Disclination model of high angle grain boundaries, Surf. Sci. 31 (1972) 12–26, http://dx.doi.org/10.1016/0039-6028(72)90251-8.

[74] A.E. Romanov, V. Vladimirov, Disclinations in crystalline solids, in: Dislocations in Solids (Vol. 9), Elsevier Science Publishers B.V., 1992, pp. 191–402, (Chapter 47).

[75] J.A. Hurtado, B.R. Elliott, H.M. Shodja, D.V. Gorelikov, C.E. Campbell, H.E. Lippard, T.C. Isabell, J. Weertman, Disclination grain boundary model with plastic dislocations by dislocations, Mater. Sci. Eng.: A 190 (1) (1995) 1–7, http://dx.doi.org/10.1016/0921-5093(94)09592-K.

[76] M. Seefeldt, L. Delannay, B. Peeters, E. Aernoudt, P. Van Houtte, Modeling the initial stage of grain subdivision with the help of a coupled substructure and texture evolution algorithm, Acta Mater. 49 (2001) 2129–2143, http://dx.doi.org/10.1016/S1359-6454(01)00126-4.

[77] A.E. Romanov, A.L. Kolesnikova, Application of disclination concept to solid structures, Prog. Mater. Sci. 54 (6) (2009) 740–769, http://dx.doi.org/10.1016/j.pmatsci.2009.03.002.

[78] M.V. Upadhyay, L. Capolungo, V. Taupin, C. Fressengeas, Grain boundary and triple junction energies in crystalline media: a disclination based approach, Int. J. Solids Struct. 48 (22) (2011) 3176–3193, http://dx.doi.org/10.1016/j.ijsolstr.2011.07.009.

[79] C. Fressengeas, V. Taupin, L. Capolungo, Continuous modeling of the structure of symmetric tilt boundaries, Int. J. Solids Struct. 51 (2014) 1434–1441, http://dx.doi.org/10.1016/j.ijsolstr.2013.12.031.

[80] X.-Y. Sun, C. Fressengeas, V. Taupin, P. Cordier, N. Combe, Disconnections, dislocations and generalized dislocations in grain boundary ledges, Int. J. Plast. 104 (2018) 134–146, http://dx.doi.org/10.1016/j.ijplas.2018.02.003.

[81] C. Zhang, A. Acharya, S. Puri, Finite element approximation of the fields of bulk and interfacial line defects, J. Mech. Phys. Solids 114 (2018) 258–302, http://dx.doi.org/10.1016/j.jmps.2018.02.004.

[82] S. Volterra, Sur l'équilibre des corps élastiques multiplement connexes, Ann. Sci. Ecol. Norm. Sup. III 24 (1907) 401–517.

[83] K.H. Anthony, U. Essmann, A. Seeger, H. Träuble, Disclinations and the cosserat continuum with incompatible rotations, in: E. Kröner (Ed.), Mechanics of Generalized Media. Proceedings of the IUTAM Symposium on the Generalized Cosserat Continuum and the Continuum Theory of Dislocations with Applications, Berlin, 1968, pp. 355–358.

[84] K.H. Anthony, Die theorie der disklinationen, Arch. Rat. Mech. Anal. 39 (1970) 43–88, http://dx.doi.org/10.1007/BF00281418.

[85] R. deWit, Linear theory of static disclinations, in: J.A. Simmons, R. de Wit, R. Bullough (Eds.), Fundamental Aspects of Dislocation Theory, vol. 1, NBS Spec. Publ. 317, National Bureau of Standards, Washington, DC, 1970, pp. 651–680.

[86] R. deWit, Theory of disclinations: II. Continuous and discrete disclinations in anisotropic elasticity, J. Res. Nat. Bur. Stand. A Phys. Chem. 77A (1) (1973) 49–100, http://dx.doi.org/10.6028/jres.077A.003.

[87] W. Huang, T. Mura, Elastic fields and energies of a circular edge disclination and a straight screw disclination, J. Appl. Phys. 41 (1969) 5175–5179, http://dx.doi.org/10.1063/1.1658641.

[88] T. Mura, Semi-microscopic plastic distortion and disclinations, Arch. Mech. Stos. 24 (1972) 449–456.

[89] H. Schaefer, The basic affine connection in a cosserat continuum, in: E. Kröner (Ed.), Mechanics of Generalized Media. Proceedings of the IUTAM Symposium on the Generalized Cosserat Continuum and the Continuum Theory of Dislocations with Applications, Springer Verlag, Berlin, 1968, pp. 57–62.

[90] E. Kossecka, R. DeWit, Disclination dynamics, Arch. Mech. 29 (1977) 749–767.

[91] C. Fressengeas, V. Taupin, L. Capolungo, An elasto-plastic theory of dislocation and disclination fields, Int. J. Solids Struct. 48 (2011) 3499–3509, http://dx.doi.org/10.1016/j.ijsolstr.2011.09.002.

[92] M. Upadhyay, L. Capolungo, V. Taupin, C. Fressengeas, Elastic constitutive laws for incompatible crystalline media: the contributions of dislocations, G-disclinations, Phil. Mag. 93 (2013) 794–832, http://dx.doi.org/10.1080/14786435.2012.733829.

[93] V. Taupin, K. Gbemou, C. Fressengeas, L. Capolungo, Nonlocal elasticity tensors in dislocation and disclination cores, J. Mech. Phys. Solids 100 (2017) 62–84, http://dx.doi.org/10.1016/j.jmps.2017.01.003.

[94] W. Dreyer, W.H. Müller, J. Olschewski, An approximate analytical 2D-solution for the stresses and strains in eigenstrained cubic materials, Acta Mater. 136 (3–4) (1999) 171–192, http://dx.doi.org/10.1007/BF01179256.

[95] B. Jiang, The least-squares finite element method, in: Theory and Computation in Fluid Dynamics and Electromagnetics, in: Springer Series in Scientific Computation, Springer, Berlin, 1998.

[96] H. Kharouji, V. Taupin, J. Guénolé, On the atomistic origin of internal length scale in strain-gradient plasticity models: The case of grain boundary structures and energies, Acta Mater. 284 (2025) 120555, http://dx.doi.org/10.1016/j.actamat.2024.120555.

[97] E. Kröner, Kontinuumstheorie der versetzungen und eigenspannungen, in: L. Collatz, F. Loesch (Eds.), Ergebnisse der Angewandte Mathematik 5, Springer Verlag, Berlin, 1958.

[98] E. Kröner, Continuum theory of defects, in: R. Balian, et al. (Eds.), Physics of Defects, North Holland, Les Houches, Session 35, New York, 1981, pp. 215–315.

[99] R. deWit, Theory of disclinations: IV. Straight disclinations, J. Res. Natl. Bur. Stand. Sect. A: Phys. Chem. 77A (5) (1973) 607–658, http://dx.doi.org/10.6028/jres.077A.036.
[100] A. Acharya, A model of crystal plasticity based on the theory of continuously distributed dislocations, J. Mech. Phys. Solids 49 (2001) 761–785, http://dx.doi.org/10.1016/S0022-5096(00)00060-0.

[101] C. Fressengeas, Mechanics of Dislocation Fields, ISTE, Wiley, 2017.

[102] C.R. Weinberger, G.J. Tucker (Eds.), Multiscale Materials Modeling for Nanomechanics, in: Springer Series in Materials Science, vol. 245, Springer International Publishing, Cham, 2016, http://dx.doi.org/10.1007/978-3-319-33480-6.

[103] R. Hadian, B. Grabowski, J. Neugebauer, GB code: A grain boundary generation code, J. Open Source Softw. 3 (29) (2018) 900, http://dx.doi.org/10.21105/joss.00900.

[104] X.-Y. Liu, F. Ercolessi, J.B. Adams, Aluminium interatomic potential from density functional theory calculations with improved stacking fault energy, Modelling Simul. Mater. Sci. Eng. 12 (4) (2004) 665–670, http://dx.doi.org/10.1088/0965-0393/12/4/007.

[105] Z. Pan, V. Borovikov, M.I. Mendelev, F. Sansoz, Development of a semi-empirical potential for simulation of Ni solute segregation into grain boundaries in Ag, Modelling Simul. Mater. Sci. Eng. 26 (7) (2018) 075004, http://dx.doi.org/10.1088/1361-651X/aadea3.

[106] E. Bitzek, P. Koskinen, F. Gähler, M. Moseler, P. Gumbsch, Structural relaxation made simple, Phys. Rev. Lett. 97 (17) (2006) 1–4, http://dx.doi.org/10.1103/PhysRevLett.97.170201, ISBN: 0031-9007.

[107] J. Guénolé, W.G. Nöhring, A. Vaid, F. Houllé, Z. Xie, A. Prakash, E. Bitzek, Assessment and optimization of the fast inertial relaxation engine (FIRE) for energy minimization in atomistic simulations and its implementation in LAMMPS, Comput. Mater. Sci. 175 (2020) 109584, http://dx.doi.org/10.1016/j.commatsci.2020.109584, arXiv:1908.02038 Publisher: Elsevier.

[108] R. Hill, The elastic behaviour of a crystalline aggregate, Proc. Phys. Soc. Sect. A 65 (5) (1952) 349–354, http://dx.doi.org/10.1088/0370-1298/65/5/307.

[109] R. deWit, Elastic constants and thermal expansion averages of a nontextured polycrystal, J. Mech. Mater. Struct. 3 (2) (2008) 195–212, http://dx.doi.org/10.2140/jomms.2008.3.195.

[110] A. Torres-Sánchez, J.M. Vanegas, M. Arroyo, Examining the mechanical equilibrium of microscopic stresses in molecular simulations, Phys. Rev. Lett. 114 (25) (2015) 258102.

[111] V. Taupin, L. Capolungo, C. Fressengeas, Disclination mediated plasticity in shear-coupled boundary migration, Int. J. Plast. 53 (2014) 179–192, http://dx.doi.org/10.1016/j.ijplas.2013.08.002.

[112] M. Frigo, S.G. Johnson, An adaptive software architecture for the FFT, Proc. Int. Conf. Acoust. Speech Signal Process. 3 (1998) 1381–1384, http://dx.doi.org/10.1109/ICASSP.1998.681704.

[113] F. Willot, Fourier-based schemes for computing the mechanical response of composites with accurate local fields, C. R. Mec. 343 (2015) 232–245, http://dx.doi.org/10.1016/j.crme.2014.12.005.