...

Wen2012.pdf

by user

on
Category: Documents
9

views

Report

Comments

Transcript

Wen2012.pdf
An image-based method for modeling the elasto-plastic
behavior of polycrystalline microstructures based on the
fast Fourier transform
Bin Wen, Nicholas Zabaras∗
Materials Process Design and Control Laboratory, Sibley School of Mechanical and
Aerospace Engineering, 101 Frank H.T. Rhodes Hall, Cornell University, Ithaca, NY
14853-3801, USA
Abstract
An efficient full-field method of computing the local and homogenized
macroscopic responses of elasto-plastic polycrystalline microstructures based
on the fast Fourier transform (FFT) algorithm is presented. This approach
takes realistic microstructure images as the input and estimates the mechanical response/properties of polycrystal microstructures under periodic boundary conditions without requiring complex discretization. Effective stressstrain response, local mechanical response fields and crystallographic texture
of deformed microstructures are examined. Special interest is given to the
fatigue indicator parameters (FIPs) of IN100 (nickel-based superalloy). This
approach accounts for both intergranular and intragranular interactions of an
elasto-plastic heterogeneous medium and therefore provides accurate prediction of the mechanical response and texture evolution. The elastic and plastic
responses are computed separately by satisfying two forms of the equilibrium
equations. A multi-grid strategy is also adopted to capture the heterogeneous
deformation of microstructures. The obtained results are compared with a
widely used crystal plasticity finite element approach. The gained computational efficiency of full-field polycrystalline microstructure simulation is
prominent.
Keywords:
Polycrystalline microstructure, Full-field simulation, Fast Fourier transform,
∗
Corresponding author: Fax:
http://mpdc.mae.cornell.edu/
607-255-1222, Email:
Preprint submitted to International Journal of Plasticity
[email protected], URL:
March 16, 2012
Form Approved
OMB No. 0704-0188
Report Documentation Page
Public reporting burden for the collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and
maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information,
including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington
VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to a penalty for failing to comply with a collection of information if it
does not display a currently valid OMB control number.
1. REPORT DATE
3. DATES COVERED
2. REPORT TYPE
16 MAR 2012
00-00-2012 to 00-00-2012
4. TITLE AND SUBTITLE
5a. CONTRACT NUMBER
An image-based method for modeling the elasto-plastic behavior of
polycrystalline microstructures based on the fast Fourier transform
5b. GRANT NUMBER
5c. PROGRAM ELEMENT NUMBER
6. AUTHOR(S)
5d. PROJECT NUMBER
5e. TASK NUMBER
5f. WORK UNIT NUMBER
7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)
Cornell University,Sibley School of Mechanical and Aerospace
Engineering,Ithaca,NY,14853
9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)
8. PERFORMING ORGANIZATION
REPORT NUMBER
10. SPONSOR/MONITOR’S ACRONYM(S)
11. SPONSOR/MONITOR’S REPORT
NUMBER(S)
12. DISTRIBUTION/AVAILABILITY STATEMENT
Approved for public release; distribution unlimited
13. SUPPLEMENTARY NOTES
14. ABSTRACT
An efficient full-field method of computing the local and homogenized macroscopic responses of
elasto-plastic polycrystalline microstructures based on the fast Fourier transform (FFT) algorithm is
presented. This approach takes realistic microstructure images as the input and estimates the mechanical
response/properties of polycrystal microstructures under periodic boundary conditions without requiring
complex discretization. Effective stressstrain response, local mechanical response fields and
crystallographic texture of deformed microstructures are examined. Special interest is given to the fatigue
indicator parameters (FIPs) of IN100 (nickel-based superalloy). This approach accounts for both
intergranular and intragranular interactions of an elasto-plastic heterogeneous medium and therefore
provides accurate prediction of the mechanical response and texture evolution. The elastic and plastic
responses are computed separately by satisfying two forms of the equilibrium equations. A multi-grid
strategy is also adopted to capture the heterogeneous deformation of microstructures. The obtained results
are compared with a widely used crystal plasticity finite element approach. The gained computational
efficiency of full-field polycrystalline microstructure simulation is prominent.
15. SUBJECT TERMS
16. SECURITY CLASSIFICATION OF:
a. REPORT
b. ABSTRACT
c. THIS PAGE
unclassified
unclassified
unclassified
17. LIMITATION OF
ABSTRACT
18. NUMBER
OF PAGES
Same as
Report (SAR)
54
19a. NAME OF
RESPONSIBLE PERSON
Standard Form 298 (Rev. 8-98)
Prescribed by ANSI Std Z39-18
Elasto-plastic response, Texture evolution, Fatigue properties
1. Introduction
Numerical prediction of effective and local mechanical behavior of polycrystalline materials has received great attention in the computational materials science community starting with the pioneering work on crystal plasticity using physically-based models by Sachs [1] and Taylor [2]. The Sachs
model, assumes that each grain in the polycrystalline aggregate is subjected
to the same stress and provides lower bound prediction to the effective stressstrain response. This model was extended to deal with large visco-plastic and
elasto-viscoplastic deformations and is known as the lower bound model [3].
On the other hand, the Taylor model also known as the upper bound model
assumes uniform strain in the polycrystal (equal to the macroscopic value),
and provides rigid prediction to the effective stress-strain response of the
microstructure. Extensions of the Taylor model to elasto-plastic [4], viscoplastic [5], and finite elasto-viscoplastic deformations [6] were proposed. Although not accounting for the interactions inside the microstructure, the
Taylor model has been widely used in computational crystal plasticity community for its simplicity and high computational efficiency [7, 8, 9, 10].
The above models, however, only give approximate (or more precisely,
extreme) estimation of the effective mechanical behavior of polycrystalline
microstructures by assuming homogeneous fields over the domain. The local
mechanical response of deformed microstructures cannot be accurately captured. To account for intergranular interaction during deformation, a more
sophisticated and popular homogenization strategy, self-consistent method,
has been developed. The formulation of self-consistent models is based on
the solution of the problem of an ellipsoidal inclusion embedded in an infinite
homogeneous equivalent medium. The inclusion is taken to be an individual
grain while the homogeneous medium represents the equivalent polycrystalline aggregate. Each inclusion (or grain) is taken as an averaged medium,
the heterogeneity within which is not considered. The first attempt to model
the overall elasto-plastic behavior of polycrystals through the self-consistent
approach was proposed by Kröner [11] based on the use of Eshelby’s solution [12]. This model neglects plastic interactions between the inclusion and
the surrounding matrix. An improvement was introduced by Hill [13] to account for the plastic interaction using an incremental formulation based on
2
the linearization of the local constitutive equations. This model was later implemented for polycrystals [14]. Hill’s incremental approach, which is known
as the secant model, was extended to study polycrystals at large elasticplastic deformations in [15]. A non-incremental self-consistent approach by
using the tangent modulus-based formulation was proposed in [16]. This work
has been used to simulate large visco-plastic deformations in various crystalline materials with full anisotropic overall tangent modulus [17, 18]. These
works were further extended to a general elasto-viscoplastic model in [19]. An
“affine” model, which yields softer predictions, was later proposed [20, 21]
to study rate-dependent elasto-plastic behavior of polycrystals. All these
developments are based on linearization schemes that, at grain level, only
make use of information on field averages, disregarding higher-order statistical information inside the grain. Improvements taking consideration of the
second-order moment of the field fluctuation in grains were also developed
in [22, 23, 24] to account for intragrain heterogeneities.
Self-consistent models are widely adopted to study the effective mechanical response and texture evolution of polycrystalline microstructures. Although they give more reliable predictions to macroscopic properties than
the Sachs and Taylor models, these “mean-field” methods are still not capable of estimating local micromechanical fields appropriately since a homogeneous grain is assumed. To this end, full-field simulations that interrogate
polycrystalline microstructures with intracrystalline resolution are developed
with the evolution of computing power. The most popular full-field model
for polycrystals nowadays is the crystal plasticity finite element method, in
which a variational solution is achieved for the force equilibrium and displacement compatibility using the principle of virtual work for a domain that is
discretized into finite elements. Several crystal plasticity finite element models have been developed incorporating different physics-based constitutive
laws since the early work by Peirce et al. [25]. For a detailed overview, the
reader is referred to [26], where theories and applications of crystal plasticity
finite element models are reviewed. In our earlier work, a rate-independent
constitutive model [27] along with the grain size effect described by lattice
incompatibility [28] were implemented to study the mechanical response of
polycrystalline FCC microstructures discretized by conforming grids [29]. A
rate-dependent nickel-based superalloy constitutive model [30] was later embedded in our crystal plasticity finite element framework to study fatigue
properties and their variability induced by microstructure uncertainties [31].
The difficulty in finite element meshing coupled with the large number of
3
degrees of freedom, however, limits the size of the microstructures that can be
investigated by finite element simulations. In addition, the high computation
cost of crystal plasticity finite element analysis also prevents it from being
used as the deterministic solver in stochastic simulations or as a point simulator in multiscale simulations. Another efficient full-field approach based on
Green’s functions and fast Fourier transform (FFT) has been proposed as an
alternative to the finite element method for solving the governing equations
for periodic heterogeneous media [32, 33]. This fast Fourier transform-based
method was originally proposed as a fast algorithm to compute the elastic
and elasto-plastic effective and local response of composites with isotropic
components [34, 35, 36]. This approach was also adapted to deal with polycrystalline microstructures using a visco-plastic constitutive model [32]. Local and macro mechanical responses, as well as the texture evolution, of 2D
and 3D realistic polycrystals were studied in a series of works [32, 37, 38].
Comparison with the self-consistent method [39] and finite element simulations [40, 41] were conducted and showed the accuracy and efficiency of the
fast Fourier transform scheme. Recent attempts to couple the fast Fourier
transform-based model with finite elements is presented in [42]. The major
merit of the crystal plasticity fast Fourier transform-based method lies in the
following aspects:
1. It directly works on pixelized microstructure images without requiring
sophisticated discretization.
2. It is a full-field method that accurately investigates the global and
local mechanical behavior of microstructures by accounting for both
intergrain and intragrain interactions.
3. It has better numerical performance than the crystal plasticity finite
element method for the same spatial resolution without sacrificing accuracy [40].
Recent crystal plasticity fast Fourier transform models neglect elastic
behavior and consider only visco-plastic constitutive behavior. For a microstructure under large plastic deformation, this approximation is valid as
elastic strain is usually very small. In certain cases, however, both elastic
and plastic mechanisms are important, e.g. fatigue analysis of turbine engine
components under cyclic loading. Driven by this requirement, we here introduce a new fast Fourier transform-based crystal plasticity model with the incorporation of elasto-plastic constitutive relations (we call it, in abbreviation,
CEPFFT). The total strain rate is additively decomposed into an elastic part
4
and a plastic part. The elastic and plastic responses are computed separately
using the fast Fourier transform-based method and combined to update the
stress field. Through a series of benchmark examples, it is shown that the
effective and local mechanical responses predicted by CEPFFT are in good
agreement with the crystal plasticity finite element results. A constitutive
model of IN100 [30] is also employed to study the microstructure-sensitive
fatigue indicator parameters of such Ni-based superalloy under cyclic loading. Comparison with finite element solutions shows great consistence of
the two methods. Besides, we analyzed the performance of the fast Fourier
transform based simulator with pure visco-plastic model implemented in two
different ways (basic formulation and augmented Lagrangian approach). It is
observed that both algorithms work equally well for crystal plasticity problems. A multi-grid strategy separating the computation and material grids
based on the particle-in-cell method [43] is also employed and the obtained
results are compared with those using a single grid strategy [40].
The organization of the paper is as follows. In Section 2, the formulation
of CEPFFT algorithm is introduced. We start with pure elastic and viscoplastic situations and extend to elasto-plastic problem. A complete boundary
value problem along with the solution strategy and elasto-plastic constitutive
relation is described. Specific constitutive models that are adopted in the
paper are presented in Section 3. The input model and microstructure update
strategy are described in Section 4. Facts about the crystal plasticity finite
element method are briefly reviewed in Section 5. Numerical examples are
demonstrated in Section 6, where comparisons between the various methods
for several benchmark problems are shown. Error analysis and computational
efficiency tests are conducted. The multi-grid and single-grid results are also
presented. Finally, a brief summary is provided in Section 7.
2. Crystal Elasto-Plastic Fast Fourier Transform Simulator
In this section, we address the solution of the boundary value problem
defining the deformation of elasto-plastic polycrystalline microstructures using Green’s function method, in which any point in the domain can be considered as an inclusion embedded in a homogeneous reference medium. The
local mechanical response of the heterogeneous medium can be calculated as
a convolution integral between the Green’s function associated with the linear reference homogeneous medium and the actual heterogeneous field [33].
All the local quantities can be written as the summation of a mean value
5
and a fluctuation indicated by a “˜” symbol. Usually, the representative volume element of bulk microstructures (not on the surface) is modeled to be
periodic and periodic boundary conditions are applied to control its deformation. In this case, the Fourier transform can be employed to efficiently
solve the problem in Fourier space, where the convolution in the real space is
reduced to simple product. An iterative scheme is needed to ensure that the
solution converges to the micromechanical responses satisfying equilibrium
and compatibility conditions.
For clarity, we will start with brief presentations of the solutions of pure
elastic and pure visco-plastic problems. Then, the solution strategy of elastoplastic problems will be introduced. For elasto-plastic problems, the strain
rate is coupled with stress and stress rate leading to challenges in using the
Green’s function method. To address this, we will introduce a solution strategy in which the elastic and plastic responses are computed simultaneously
but separately.
2.1. Solution of Crystal Elastic Boundary Value Problems
In pure crystal elastic problems, the total strain rate ε̇ is equal to the
elastic strain rate ε̇e and the stress rate is linearly related to the strain rate
through the generalized Hooke’s law:1
σ̇(x) = Ce (x) : ε̇e (x),
(1)
e
e
is the local elastic stiffness
where Cijkl
(x) = Rim (x)Rjn (x)Rko (x)Rlp (x)Čmnop
tensor represented in the sample coordinate system that relates the crystal
lattice frame via a rotation matrix R(x) determined by the orientation r(x)
of the crystal at position x. Če is the elastic stiffness modulus in the lattice
coordinate system.
For a microstructure subjected to periodic boundary conditions with an
imposed average velocity gradient L = ∇V, the local equilibrium equation,
represented in terms of stress rate, needs to be satisfied at any point x within
the microstructure domain B. The complete boundary value problem is defined as:
ṽe
˜
∇ · σ̇(x) = ∇ · σ̇(x)
=0
∀x ∈ B,
is periodic,
σ̇ · n is antiperiodic on ∂B,
1
(2)
We state the elastic problem in terms of stress and strain rates for exploring similarities
with the formulation of the elasto-plastic problem in Section 2.2.
6
where ṽe (x) = ve (x) − L · x is the velocity fluctuation (deviation of the local
velocity ve (x) from the mean velocity V) at x induced by the microstructure
heterogeneity. The superscript e indicates that the response stems from
elastic deformation. Our goal is to compute the velocities and their gradients
of all material points that satisfy the above governing equations and use them
to evaluate the strain and stress responses over the microstructure domain.
To this end, we can write the local stress σ(x) as the sum of two terms:
σ̇(x) = C0 : ε̇e (ve (x)) + φe (x).
(3)
In Eq. (3), C0 is the stiffness modulus of a linear homogeneous reference
medium, in which point x is imbedded. In this elastic problem, C0 is selected
as the averaged elastic modulus Ce over the microstructure domain:
Z
1
0
e0
e
C = C = hC ih =
Ce (x)dx.
(4)
VB B
The second term, φe (x), in Eq. (3) is the periodic polarization field defined as
˜
˜e (ṽe (x)).
φe (x) = σ̇(x) − Ce0 : ε̇e (ve (x)) = σ̇(x)
− Ce0 : ε̇
(5)
Substituting Eq. (3) into the equilibrium equation (Eq. (2)), we obtain
e0 e
Cijkl
ε̇kl,j + φeij,j = 0
e0 ˜e
Cijkl
ε̇kl,j + φeij,j = 0.
or
Representing the strain rate in terms of the velocity gradient ε̇eij =
we obtain:
e0 e
Cijkl
vk,lj + φeij,j = 0
e0 e
Cijkl
ṽk,lj + φeij,j = 0.
or
(6)
1
2
e
e
vi,j
+ vj,i
,
(7)
The differential Eq. (7) can be solved by means of the Green’s function
method. Introducing the Green’s function Ge (x, x′ ), the solution ṽke (x) takes
the form
Z
e
ṽk (x) = − Gekm (x − x′ )φemn,n (x′ )dx′ .
(8)
B
Substituting Eq. (8) into Eq. (7), leads to:
Z
Z
e0
e
′
e
′
′
−Cijkl Gkm,lj (x − x )φmn,n (x )dx + δim φemj,j (x′ )δ(x − x′ )dx′ = 0, (9)
B
B
7
which can be rearranged as
Z
e0
−Cijkl
Gekm,lj (x − x′ ) + δim δ(x − x′ ) φemn,n (x′ )dx′ = 0.
(10)
B
Taking φemn,n to be arbitrary, we arrive at the local equilibrium equation in
the Green’s function form:
e0
−Cijkl
Gekm,lj (x − x′ ) + δim δ(x − x′ ) = 0.
(11)
This equation can be transformed to Fourier space where the convolutional
solution in the real space (Eq. (8)) is represented by a simple product. The
equilibrium equation in Fourier space is then:
e0
ξl ξj Cijkl
Ĝekm (ξ) = −δim ,
(12)
where ξ is a point (frequency) in the Fourier space. Solving the linear system
Eq. (12), we obtain the Green’s function in Fourier space Ĝekm to be
e0
Ĝe = Ae−1, Aeik = −ξl ξj Cijkl
.
(13)
Γ̂eijkl = −ξl ξj Ĝeik ,
(14)
Defining
and integrating Eq. (8) by parts while assuming that the boundary terms
vanish [32], we can compute the velocity fluctuation as:
Z
e
ṽk (x) =
Gekm,n (x − x′ )φemn (x′ )dx′ .
(15)
B
The velocity and its gradient fluctuations2 in Fourier space are
e
ṽˆi (ξ) = iξj Ĝeim (ξ)φ̂emj (ξ),
e
ṽˆi,k (ξ) = Γ̂eikmj (ξ)φ̂emj (ξ).
2
(16)
(17)
Note that the fluctuations can also be total velocity and its gradient, depending on
the assignment of the 0 frequency term.
8
ˆ e (ξ))),
After transforming them back to the real space (e.g. ṽe (x) = F F T −1 (ṽ
the strain rate and spin rate fluctuations can be evaluated, respectively, by
e
e
˜ε̇eij (x) = 1 (ṽi,j
+ ṽj,i
),
2
e
e
˜ eij (x) = 1 (ṽi,j
− ṽj,i
).
ω̇
2
(18)
The stress rate can be updated according to the Hooke’s Law (Eq.(1)). With
the updated stress rate and strain rate, we can perform the next iteration
until converged results are reached.
2.2. Solution of Crystal Visco-Plastic Boundary Value Problems
If the deformation is assumed to be pure visco-plastic (i.e. the elastic
response is completely neglected, ε̇ = ε̇p ). A nonlinear constitutive model is
employed to link stress to strain rate in the form:
ε̇p (x) = Mp (σ(x)) : σ(x),
(19)
where Mp (σ(x)) is a plastic compliance tensor that is nonlinearly dependent
on stress σ(x).
The solution procedure of a visco-plastic problem is similar to that in
Section 2.1 except that this time the equilibrium equation is written in terms
of stress rather than stress rate:
ṽ
p
∇ · σ(x) = ∇ · σ̃(x) = 0
∀x ∈ B,
is periodic,
σ · n is antiperiodic on ∂B,
(20)
where the superscript p indicates plastic deformation induced response. Incompressibility of plastic deformation also needs to be satisfied. In [32], a
Fourier transform-based algorithm was proposed where the incompressibilp
ity condition was satisfied by introducing explicitly the constraint vk,k
= 0.
In the current work, the incompressibility condition is accounted for by
computing the polarization φp (x) with a strain rate updated iteratively as
ε̇p ← ε̇p − 13 tr(ε̇p ) as follows:
˜p (ṽp (x)),
φp (x) = σ(x) − Cp0 : ε̇p (vp (x)) = σ̃(x) − Cp0 : ε̇
(21)
where the stiffness modulus of the linear homogeneous reference medium is
taken to be the averaged plastic modulus Cp = Mp−1 over the microstructure
9
domain:
C
p0
1
= hC ih =
VB
p
Z
Cp (σ(x))dx.
(22)
B
The form of the plastic modulus Cp or equivalently the plastic compliance Mp
is determined by the specific plastic constitutive model that is adopted. In
Section 3, we will present two constitutive models to be used in the examples.
The plastic problem can be solved following the steps from Eq. (6) to Eq. (18)
simply by replacing the superscript e to p and using stress instead of stress
rate.
2.3. Solution of Crystal Elasto-Plastic Boundary Value Problems
For elasto-plastic problems, the total strain rate ε̇ is additively decomposed into an elastic term ε̇e and a plastic part ε̇p with tr(ε̇p ) = 0:
ε̇(x) = ε̇e (x) + ε̇p (x).
(23)
Following the constitutive relations given above, both stress σ and stress rate
σ̇ are entangled with strain rate:
ε̇(x) = Me (x) : σ̇(x) + Mp (σ(x)) : σ(x),
(24)
where Me = Ce−1 is the elastic compliance tensor. The local mechanical
response at time t depends on the entire loading history of the specimen.
The current stress σ is computed by σ = σ n + σ̇dt, where σ n is the stress
at the previous time step.
To follow a Green’s function approach to the elasto-plastic boundary value
problem, a proper modulus C0 , taking both elasticity and plasticity into account, needs to be designed for the homogeneous reference medium that can
directly link either stress or stress rate to strain rate. However, to design such
an effective modulus is not trivial. Therefore, we here propose a formulation
that solves for the elastic and plastic velocity fluctuations separately using
the fast Fourier transform-based algorithm thus avoiding the construction of
the elasto-plastic modulus. The total velocity gradient at a single point is
then computed by adding the two fluctuations to the mean value ∇V. After
that, a nonlinear constitutive model is designed to update the stress and
stress rate given the total strain rate. We will this approach as CEPFFT.
The key of the CEPFFT approach is that we solve simultaneously the two
forms of the equilibrium equations defined in Eq. (2) and (20) for elastic and
10
plastic velocity (and their gradient) fluctuations. The total velocity gradient
can then be obtained as
∇v(x) = ∇V + ∇ṽe (x) + ∇ṽp (x),
(25)
from which, the total strain rate ε̇ can be calculated as the symmetric part
of ∇v. However, the portion of elastic strain rate ε̇e and plastic strain rate
ε̇p in ε̇ is not known. As a result, the stress and stress rate corresponding to
a given total strain rate cannot be directly computed.
An iterative scheme is designed to linearize the nonlinear relations among
stress, stress rate, and strain rates, in order to provide a way of finding the
elastic and plastic strain rates given the total strain rate. We first rewrite
Eq. (23) as
F = ε̇e (x) + ε̇p (x) − ε̇(x) = 0.
(26)
We aim at solving this equation for the elastic strain rate ε̇e with known ε̇
using the Newton-Raphson scheme. To this end, Eq. (26) can be linearized
as follows:
dF
(27)
F(i+1) (ε̇e(i+1) ) = F(i) (ε̇e(i) ) + e ε̇e(i+1) − ε̇e(i) .
dε̇
According to the elastic and plastic constitutive models, as well as the stress
incremental equation σ = σ n + ∆tσ̇, we can write the following relations:
dε̇p
= Mpt ,
dσ
dσ
d (σ n + ∆tσ̇)
=
= ∆tII,
dσ̇
dσ̇
dσ̇
e
e = C ,
dε̇
(28)
where II is the fourth-order identity tensor and the tangent plastic compliance Mpt is specified by the particular plastic constitutive model used (see
Section 3.1). Using these relations, we can simplify Eq. (27) as:
e
dε̇p
dε̇
e(i+1)
e(i)
(i+1) e(i+1)
(i) e(i)
:
ε̇
−
ε̇
+
F
(ε̇
) = F (ε̇ ) +
dε̇e dε̇e
e
dε̇p dσ dσ̇
dε̇
(i) e(i)
: ε̇e(i+1) − ε̇e(i)
= F (ε̇ ) +
e +
e
dε̇
dσ dσ̇ dε̇
(i) e(i)
= F (ε̇ ) + (II + ∆tMpt : Ce ) : ε̇e(i+1) − ε̇e(i) . (29)
11
Setting F(i+1) = 0, the elastic strain rate at iteration i + 1 is computed
by
ε̇e(i+1) = ε̇e(i) − (II + ∆tMpt : Ce )−1 : F(i) .
(30)
With the elastic strain rate ε̇e(i+1) at the (i+1)th iteration computed from
the equation above, the stress rate σ̇ (i+1) can be obtained using Eq. (1).
Therefore, the stress is updated as σ (i+1) = σ n + σ̇ (i+1) , with which the
plastic strain rate ε̇p(i+1) can be computed using the plastic constitutive relation Eq. (19). The impressibility is enforced by setting ε̇p(i+1) ← ε̇p(i+1) −
1
tr(ε̇p(i+1) ). Iteratively updating the above equations, the final elastic as
3
well as the plastic strain rates can be computed until convergence is achieved
(i.e. when F in Eq. (26) approaches 0). The converged stress rate and stress
can be subsequently computed according to the constitutive model. After
that, we can construct the elastic and plastic polarization fields, respectively,
following Eqs. (5) and (21). By transforming them to Fourier space, the
fluctuations of velocity gradients induced by elastic and plastic deformation
can be updated using Green’s functions. Inversely transforming these fluctuations to real space, a new strain rate field as well as stress and stress rate
can be obtained. The algorithm can then proceed to the next iteration. The
overall CEPFFT algorithm is summarized next.
2.4. CEPFFT Algorithm
We adopt a basic fast Fourier transform-based algorithm for implementing the CEPFFT simulator. The algorithm was firstly proposed in [34, 35]
to deal with elastic and elasto-plastic heterogeneous composites with moderate contrast. It was also adapted by Lebensohn [32] to solve elastic and
visco-plastic polycrystalline problems. This method is based on the exact
expression of Green’s function for linear elastic, homogeneous reference material.
Algorithm:
1. At the first iteration, start with an initial guess of the total velocity
gradient field at time step n + 1: 0 ∇vn+1 (x) = ∇vn (x), ∀x ∈ B, from
which the local strain rate can be computed ε̇(x) = sym(∇v(x)). Then
evaluate the elastic portion ε̇e (x) and plastic portion ε̇p (x) of the strain
rate using the local constitutive relations. At the same time, compute
the initial stress 0 σ(x) and stress rate 0 σ̇(x).
12
2. Compute the elastic and plastic polarization fields, i φe (x) and i φp (x),
respectively, for the i iteration given the stress, stress rate and strain
rate fields:
φe (x) =i σ̇(x) − Ce0 :i ε̇e (x),
i p
φ (x) =i σ(x) − Cp0 :i ε̇p (x).
i
(31)
3. Transform the polarizations to Fourier space via fast Fourier transform:
p
i e
φ̂ (ξ) = F F T (i φe (x)) and i φ̂ (ξ) = F F T (i φp (x)).
4. Compute the velocity gradient fluctuation fields induced by elastic and
plastic deformation, respectively, in the Fourier space for the (i + 1)-th
iteration.
e
e
e
ˆ (ξ) = Γ̂ (ξ) :i φ̂ (ξ), ∀ξ 6= 0; and
∇ṽ
p
p
i+1 ˆ p
∇ṽ (ξ) = Γ̂ (ξ) :i φ̂ (ξ), ∀ξ 6= 0; and
i+1
e
ˆ (0) = 0,
∇ṽ
i+1 ˆ p
∇ṽ (0) = 0.
i+1
(32)
5. Transform the current velocity gradient fluctuation fields back to the
e
real space
i.e. i+1
through
inverse fast Fourier transform,
∇ṽ (x) =
e
p
ˆ (ξ) and i+1 ∇ṽp (x) = F F T −1 i+1 ∇ṽ
ˆ (ξ) .
F F T −1 i+1 ∇ṽ
˜e +i+1 ε̇
˜p and then the
6. Compute the total strain rate i+1 ε̇ = Ė +i+1 ε̇
stress i+1 σ(x) and stress rate i+1 σ̇(x) fields according to the constitutive model.
7. Check the error (equilibrium condition):
1/2
hk∇ · i+1 σk2 i
e=
ki+1 σk
hkξ · i+1 σ̂k2 i
=
ki+1 σ̂(0)k
1/2
.
(33)
If e is smaller than a prescribed error tolerance, the iteration process
stops. Otherwise, return to step (2) and proceed to the next iteration.
Upon convergence, the grain orientations and the positions of the pixel
points are updated according to the spin and velocity gradient fields, respectively (see Section 4.3 for details).
Remark 1: The error of the equilibrium condition (Step (7)) is mostly
determined by the resolution of the image (as will be shown later in the
examples). For images with coarse resolution, the error may converge, with
fluctuations, to a value larger than 0. This error is inherently associated
13
with the FFT-based methodology. Therefore, a practical way to check the
convergence is to examine the relevant difference between the error at the
current and last iterations:
η=
|i+1 e −i e|
.
i+1 e
(34)
Remark 2: It is clear that the computational cost of one iteration step of the
algorithm discussed above is doubled that corresponding to either the pure
elastic or visco-plastic algorithms. If a proper elasto-plastic modulus can be
designed, the problem can be solved using only one set of equations, which
may lead to higher computational efficiency. A feasible approach arises by
representing the local stress σ(x) in terms of the total strain rate ε̇(v(x)):
σ(x) = C0 : ε̇(v(x)) + φ(x),
(35)
where C0 = hCep ih is computed as the volume average of an elasto-plastic
modulus Cep derived as follows:
Cep =
dσ dε̇e
dσ dε̇p
dσ
= e
+ p
= ∆tCe + Cp .
dε̇
dε̇ dε̇
dε̇ dε̇
(36)
The plastic modulus Cp is chosen to be the secant modulus that depends
on the specifics of the constitutive model adopted. With the construction
of the elasto-plastic modulus, the problem can be solved following the same
procedure as for the visco-plastic problem. The convergence rate of this
particular algorithm will be shown in Section 6 that is slower than the main
algorithm presented earlier. A more sophisticated design of the elasto-plastic
modulus is of great interest.
3. Constitutive Model
The fast Fourier transform based full field algorithm can be adapted to
various materials assuming appropriate constitutive models are provided. In
this section, we will describe two plastic constitutive models that are adopted
in the numerical examples for single phase FCC crystal and IN100 superalloy.
The plastic modulus that is needed in the CEPFFT method is derived. The
fatigue indicator parameters that will be employed to measure the fatigue
properties of Ni-based superalloys are defined in [44].
14
3.1. Plastic Constitutive Relations
Various plastic constitutive relations have been developed based on different flow rules and hardening laws. Usually, the nonlinear connection between
plastic strain rate and Cauchy stress, ε̇p (σ(x)), is assumed, where the plastic strain rate ε̇p is defined as the sum of the contributions of shear rates,
γ̇ (α) (x), on all active slip systems:
p
ε̇ (x) =
Ns
X
m(α) (x)γ̇ (α) (x).
(37)
α
In the above equation, Ns is the number of active slip systems, and m(α)
denotes the symmetric Schmid tensor of slip system α:
m(α) = sym(S(α) ) =
1 (α)
s ⊗ n(α) + n(α) ⊗ s(α) ,
2
(38)
where s(α) and n(α) are the slip direction and slip plane normal of the system
α, respectively.
The shear rate γ̇ (α) (x) is determined by the resolved shear stress τ (α) =
S(α) : σ = m(α) : σ and slip resistance κ(α) on slip system α through a specific
flow rule. The difference between the various plastic constitutive models lies
in the choice of the flow rule and hardening law that defines the evolution of
the slip resistances.
3.1.1. A typical flow rule and hardening law for single phase crystalline materials
The first model we introduce here is a typical one for single phase polycrystals. Asaro and Needleman [6] proposed a widely used rate-dependent
crystal plasticity flow rule for evaluating shearing rate, γ̇ (α) , on slip system
α:
(α) (1/m)
τ (α)
γ̇ = γ̇0 (α) sign(τ (α) ),
(39)
κ
where γ̇0 is a reference rate of shearing and m characterizes the material rate
sensitivity. κ(α) is the slip resistance of system α. The slip resistance evolves
following certain hardening law that is calibrated from experimental data to
15
describe the strain hardening effect. Currently, we adopt a Voce type model
given in [40], where the hardening rate κ̇(α) is computed by
κ̇(α) =
dκ̄(α) X αβ (β)
h γ̇ ,
dΓ β
(40)
with the hardening function
(α)
κ̄
θ0 Γ
= κ0 + (κ1 + θ1 Γ) 1 − exp −
,
κ1
(41)
where hαβ is a hardening matrix whose diagonal elements denote self-hardening
and off-diagonal elements denote latent hardening. In the current work, we
assume both of the two are identical and equal to 1. κ0 , κ1 , θ0 , θ1 are material
dependent parameters. The cumulative shear Γ in a grain is defined as
Z tX
Γ=
γ̇ (α) ∆t.
(42)
0
α
Following the above equations, we can derive an explicit form of the
plastic constitutive law as
ε̇p (x) = Mps (σ(x)) : σ(x),
(43)
where Mps is called the secant plastic compliance having the form given below:
(1/m−1)
Ns
X
m(α) (x) ⊗ m(α) (x) m(α) (x) : σ(x) p
Ms (σ(x)) = γ̇0
.
(44)
κ(α) (x)
κ(α) (x)
α
The corresponding tangent plastic compliance in Eq. (30) at this specific case is Mpt = m1 Mps . The plastic modulus in Eq. (36) for the method
highlighted in Remark 2 is Cp = Mp−1
s . Given a plastic strain rate, the
corresponding stress needs to be computed following a iterative scheme (e.g.
Newton-Raphson method) because of the nonlinear nature of the constitutive
model.
3.1.2. A specific flow rule and hardening law for two-phase nickel-based superalloy IN100
In this work, we are also interested in the fatigue properties of IN100, a
nickel-based superalloy, at high temperature (650◦C). Since the microstructure of superalloys is complex, more sophisticated flow rule and hardening
16
laws are developed to describe the micromechanical behavior of two-phase
microstructure. In this work, we employ the homogeneous constitutive model
developed in [30]. The second phase (γ ′ precipitates) configuration is not explicitly modeled. Instead, the effects of the second phase are taken into
account through particular parameters. Cube slip h110i{100} systems are
activated in addition to octahedral slip h110i{111} systems to take cross slip
mechanism at high temperatures into consideration. The rate dependent flow
rule which estimates the shearing rate on each slip system includes a back
force term χ for the modeling of the Baushinger effect arising principally from
matrix dislocation interaction with γ ′ phase. The effect of volume and size
of γ ′ precipitates on material strength is taken into account by constitutive
parameters. The constitutive equations are summarized below and detailed
in [30, 44, 45].
The flow rule of slip system α is
+
*
"
(α)
(α) n1
(α)
|τ
−
χ
|
−
κ
(α)
λ
λ
γ̇ (α) = γ̇1
(α)
Dλ
+n2 #
*
(α)
(α)
|τ
−
χ
|
(α)
(α)
λ
sgn(τ (α) − χλ ),
(45)
+ γ̇2
(α)
Dλ
(α)
(α)
where γ̇1 and γ̇2 are constants related to the initial shearing rate and
(α)
Dλ is the drag stress assumed to be constant. λ ={oct, cub} refers to the
octahedral and cube slip systems, respectively. The function hxi returns x if
x > 0 and returns 0, otherwise.
(α)
The evolution of the slip resistance κλ follows the Taylor strain harden(α)
ing law determined by dislocation density ρλ :
q
(α)
(α)
(α)
(46)
κλ = κ0,λ + αt µmix b ρλ ,
where αt = 0.0385 and µmix = (fp1 + fp2 + fp3 )µγ ′ + fm µγ . µγ ′ and µγ are
shear moduli for γ ′ precipitates and γ matrix, respectively. The magnitude
of Burgers vector is b = (fp1 + fp2 + fp3 )bγ ′ + fm bγ . fp1 , fp2 , fp3 are volume
fractions of primary, secondary, and tertiary γ ′ precipitates, respectively, and
p1
′
,
fm = 1−fp1 −fp2 −fp3 is the volume fraction of γ matrix phase. fp1
= fp1f+f
m
f
f
p3
p2
and fp3 = fp3 +f
. For different slip systems, the initial slip
fp2 = fp2 +f
m
m
resistance can be evaluated by
i1/nk
h
(α)
(α)
nk
(α)
+ (fp1 + fp2 )τns
,
κ0,oct = (τ0,oct )nk + ψoct
17
i1/nk
h
(α)
nk
,
(τ0,cub )nk + ψcub
(α)
κ0,cub =
where
ψλ = cp1
w =
s
′
fp1
w
+ cp2
d1
ΓAP B
ΓAP B−ref
,
s
w
(47)
q
′
fp2
cgr
′
d3 + p ,
+ cp3 wfp3
d2
dgr
(48)
and
(α)
(α)
(α)
(α)
τns
= hpe τpe
+ hcb |τcb | + hse τse
,
(49)
in which ΓAP B is the anti-phase boundary energy density here taken be equal
to ΓAP B−ref , di, i = 1, 2, 3 are the sizes of precipitates, and dgr is the grain
(α)
size. The non-Schmid stress τns is related to the resolved shear stresses on
(α)
(α)
(α)
the primary τpe , cube τcb and secondary τse slip systems, respectively, by
parameter hpe , hcb , and hse . Currently, we assume this term to be 0.
The dislocation density evolution has the following form:
q
(α)
(α)
(α)
ρ̇λ = h0 Z0 + k1,λ ρλ − k2,λ ρλ
|γ̇ (α) |,
−1
2
kδ
, dδef f ≈
.
(50)
Z0 =
bdδef f
d2δ
(α)
The evolution of the back stress χλ is also dependent on dislocation
density and shear rate:
q
(α)
(α)
(α)
(α)
(α)
|γ̇ (α) |,
(51)
χ̇λ = Cχ ηλ µmix b ρλ sgn(τ − χλ ) − χλ
ηλ =
η0,λ Z0
q
,
(α)
Z0 + k1,λ ρλ
′
′2
where Cχ = 123.93 − 433.98fp2
+ 384.06fp2
.
◦
Parameters for superalloys at 650 C listed in [30] are adopted. For additional information about the constitutive model, the reader is referred to [44].
Given the stress, the plastic strain rate needs to be computed following the
18
nonlinear relationship and vice versa. Eq. (30) requires the evaluation of
tangent plastic compliance. In the case of IN100 superalloy, the compliance
can be derived as:
Ns dε̇p dγ̇ (α) dτ (α)
dε̇p X
p
=
Mt (x) =
dσ
dγ̇ (α) dτ (α) dσ
α

*
+
Ns
(α)
(α)
(α) n1 −1
X
n1
|τλ − χλ | − κλ
=
m(α) (x) ⊗ m(α) (x) γ̇1 (α)
(α)
Dλ
Dλ
α
*
+n2 −1 
(α)
(α)
|τλ − χλ |
n2
.
(52)
+ γ̇2 (α)
(α)
Dλ
Dλ
The secant plastic modulus required by Eqs. (21) and (36) is taken to be
Cp = Mp−1
with
s

*
+
Ns
(α)
(α)
(α) n1 −1
X
1
|τ
−
χ
|
−
κ
λ
λ
λ
m(α) (x) ⊗ m(α) (x) γ̇1 (α)
Mps (x) =
(α)
Dλ
Dλ
α

*
+n2 −1
(α)
(α)
|τλ − χλ |
1
.
(53)
+ γ̇2 (α)
(α)
Dλ
Dλ
With these relations, the algorithm proposed in Section 2 can be applied to
IN100 superalloy.
4. Microstructure Model
The solution strategy of computing the deformation of polycrystalline
microstructures under periodic boundary conditions was discussed in Section 2. The main procedure is: (1) compute polarization field; (2) transform
the polarization field to Fourier space using fast Fourier transform; (3) update velocity gradient in the Fourier space and transform it back to the real
space; (4) update real-space fields (e.g. strain rate, spin rate, stress rate,
etc.) accordingly. An iterative scheme is adopted to obtain convergence.
In this section, we would introduce the digital microstructure model that is
used as the input to FFT-based simulations. The update strategy of the
microstructure and crystal orientation during deformation is also described.
19
4.1. Discretization
The input to the FFT-based (including pure elastic, pure visco-plastic,
and elasto-plastic) simulators is pixelized images with orientation parameters
associated with each pixel (or voxel for 3D). The pixels or voxels are the
discretization of the input image, which essentially requires no effort.
To implement the FFT-based algorithm and solve the underlying boundary value problem, the microstructure is discretized by a regular grid consisting of N1 × N2 pixels (2D problem) or N1 × N2 × N3 voxels (3D problem).
Denote Li to be the period (edge length) of the microstructure in the ith
direction (i = 1, 2 for 2D and i = 1, 2, 3 for 3D). The coordinates of the
point, which is pixel in 2D and voxel in 3D, in the ith direction is therefore:
xi = 0,
Li
Li Li
, 2 , . . . , (Ni − 1) .
Ni Ni
Ni
(54)
In Fig. 1, examples of 2D and 3D grids of polycrystalline microstructures
are shown in comparison with their image views. Different colors in the microstructure represent grains with distinct orientations. The grain structures
are generated following a Voronoi tessellation scheme, where the positions of
centroids are adjusted to minimize the interaction forces [46] (see Section 4.2).
The regular discretization grid of the microstructure determines a regular
reciprocal grid in the Fourier space, which makes the fast Fourier transform
convenient. The i-th direction coordinates of the points in the reciprocal
grid, namely frequencies, are
Ni
Ni
1
1
ξi = − + 1
, − +2
,...,
2
L
2
L
i
i
1
1
Ni
1
Ni 1
− , 0, , . . . ,
−1
,
,
(55)
Li
Li
2
Li
2 Li
where i = 1, 2 for 2D and i = 1, 2, 3 for 3D. In the current work, the number
of points in each dimension is selected to be a power of 2 in order to facilitate
the fast Fourier transform that is conducted using the FFTW libraries [47].
4.2. Microstructure Generation
The microstructure image can be generated by several experimental or
numerical ways, e.g. electron backscattered diffraction (EBSD), phase field
method, Monte Carlo grain growth, or Voronoi tessellation. In the current
20
Figure 1: (a) The image representation of a 2D polycrystalline microstructure containing
10 grains. (b) The pixel grid of the 10-grain 2D microstructure. The microstructure
is discretized by 16 × 16 pixels. (c) The image representation of a 3D polycrystalline
microstructure containing 64 grains. (b) The voxel grid of the 64-grain 3D microstructure.
The microstructure is discretized by 16 × 16 × 16 voxels.
work, we adopt the Voronoi tessellation scheme in combination with a “relaxation” strategy to generate polycrystalline microstructures with approximately controlled grain sizes. The algorithm is summarized as follows:
• Given a set of target sizes of the microstructure constituent grains, we
first randomly place Voronoi cell centroids within the microstructure
domain. The number of centroid points, Np , is decided by Np = V /µd ,
where V is the microstructure volume (for 3D) or area (for 2D), and
µd is the mean volume or area of grains.
• A force function F(i, j) is defined for the adjustment of grain centroid
21
locations:
F(i, j) = n(j, i) ∗ hr(i) + r(j) − dist(i, j)i,
(56)
where F(i, j)is the force exerted by point j on point i in the direction
n(j, i), which is along the line starting from j and pointing to i. r(i)
is the radius (assuming spherical (for 3D) or circular (for 2D) grains)
of grain i. dist(i, j) is the distance between point i and point j, which
can be periodic if the microstructure is periodic. Function hxi = x, if
x > 0, and hxi = 0, otherwise.
• Relaxation step: move points according to the velocity function defined
as
v(i) = F(i)/V (i),
(57)
P
where F(i) is the net force of point i: F(i) = j F(i, j) and V (i) is the
volume or area of the corresponding grain. The location of centroid i
is updated to Pn+1 (i) = Pn (i) + v(i)∆t, where ∆t adaptively changes
for the convenience of convergence.
• After all centroid locations are stable (force balance is reached), a
Voronoi tessellation is computed using them. The polycrystalline microstructure is therefore constructed.
A schematic description of relaxation of the Voronoi cell centroids is depicted in Fig. 2.
4.3. Grid and Texture Update
The grid after deformation may become irregular as material points move
according to local velocities. The exact new position of material point X is
x(X) = X + (L + ∇ṽ(X)) X∆t,
(58)
where L = ∇V is the homogeneous (average) velocity gradient of the microstructure and ṽi,j (x) is the velocity gradient fluctuation at point x.
To model the deformation of the microstructure, an irregular material grid
is expected. However, this irregular grid in the real space results in difficulties
on conducting Fourier transform in the next time step in a time-dependent
simulation, because that fast Fourier transform requires a regular grid. To
22
1
1
2 10
3
8
0.9
5
0.8
0.8
9
Relaxing
centroid
positions to
minimize the
“Energy”.
0.7
0.6
0.5
4
7
0.4
10
10 9
9 3
3 3 3
3 2 2 2 2 2 10 10 10 10
9 9
9 3
3 3 3
3 2 2
2 2
2 10 10 10
9 9
9 3
3 33 3
3 2 2
2 2
5 10 10 9
9 99
9 9
3 3 3
3 3 5
5 5
5 5 9
9
9 9
9 9
3 3 3
3 5 5
5 5
5 5 9
9
5 5 4
4
0.9
0.7
0.6
0.5
0.3
0.2
0.2
1
0.1
9 9
7 7 7
7 5 5
4 4
4 7
7 7 7
7 7 5
5 5
5 5 4
4
4 4
4 7
7 7 7
7 7 5
5 5
5 4 4
4
4 4
4 7
7
7
7 7
7 7 1
1 1
1 4 4
44
4 4
4 7
7 7 7
7 7 1
1 1
1 1 4
4
4 4
6 6
7 7 7
7 1 1
1 1
1 1 4
4
11
1 1 1
4
0.4
0.3
6
9 9
55 5
6 6
6 6
6 8 8
8 8 1
1
6 6
6 6
6 8 8
8 8 8
1 1
1 1 10 10
6 6
66
6
6 8 8
1 1
1 10 10 10
10 6
6 6
6 8 8
88
8 8
0.1
0 10
0
0
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
(a)
60.1 6 0.2
6
8 8 2
2 2 10 10 10 10
60.3 3 0.4
8 8 0.52 20.6 2 0.7
2 100.810 10
10 1
0.9
(b)
Figure 2: (a) The initial randomly placed Voronoi cell centroids with corresponding circular
grain envelopes. (b) The relaxed centroids which are used for Voronoi tessellation. The
microstructure is assumed to be periodic. Numbers in figures are grain IDs.
resolve this complexity, a strategy of using two grids, a regular computation
grid and an irregular material configuration grid, proposed in [43] based on
the Particle-In-Cell (PIC) method [48, 49] is introduced. The computation
grid is used for applying the fast Fourier transform method to evaluate the
strain related fields. It is a regular grid but not necessary rectangular as
required by Fourier transform. The material grid, on the other hand, is
attached to material particles, on which constitutive relations are carried
out.
Each grid carries its own set of unknowns. Information needs to be transferred back and forth between the two grids during computation. At the beginning of each time step, the initial guess of local stress and polarization are
computed on the material grid given the initial strain rate. The polarization
field is transferred to the computation grid, on which fast Fourier transform
is performed. The updated velocity gradient is then transferred back to the
material grid so that the stress related fields can be updated using the constitutive model. At the end of the time step, the computation and material
grids are deformed. The material particles move according to local velocity
(Eq. (58)) and the regular computation grid evolves with the average velocity
gradient:
x(X) = (I + L∆t) X.
(59)
A schematic showing the operation of the multi-grid strategy is depicted in
Fig. 3.
23
A
vi , j (x p ) = ∑a=1 a (x p )vi, j (xca )
Real space
velocity gradient
Real space
velocity gradient
invFFT
Real space
strain rate
Constitutive
Fourier space
velocity gradient
model
Green’s
Real space
stress
operator
Fourier space
polarization
FFT
Real space
polarization
Real space
polarization
Material Grid
Computation Grid
1
K
k a
φij (xc ) =
∑ mp (xkp )φij (xkp )
m(xc ) k =1
Figure 3: A schematic description of the multi-grid CEPFFT strategy. The constitutive
model is applied only on the material grid, whereas fast Fourier transform operates on the
computation grid.
The information transfer from regular computation grid to material points
is performed through an interpolation operation. A family of shape functions
N a (x), a = 1, . . . , A are introduced. Since the regular computation grid
forms quadrilateral (in 2D) or brick (in 3D) elements, a natural choice of N a
are the classical finite element shape functions. Given nodal values fc (xac )
of a function f (x) on the computation grid (denoted by subscript c), the
interpolated value at any material point xp in the microstructure (denoted
by subscript p) is given as
fp (xp ) =
A
X
N a (xp )fc (xac ).
(60)
a
The interpolation is local since finite element shape functions are adopted.
Only nodes belonging to the element that contains particle xp are involved.
The inverse transfer from the material grid to the computation grid is
performed using the same interpolation functions [43]. We assume all material particles have the same mass mp . The mass of a computational node is
24
defined as
mc (xc ) =
K
X
mkp N a (xkp ).
(61)
k
The sum is over all material particles that are contained in the elements
that share a common computation node xc . N a is the corresponding shape
function that connects particle xkp to node xc . The function value fc of f at
any computation node xc is therefore:
K
1 X k a k
mp N (xp )fp (xkp ).
fc (xc ) =
mc (xc ) k
(62)
A 2D illustration of initial and deformed computation and material grids
is shown in Fig. 4, where big red spots represent material particles and small
black spots are the nodes of the computation grid. The deformed computation grid remains regular so that fast Fourier transform can be conducted,
while material particles move heterogeneously. In most PIC studies, the material grid is finer than the computational one. In the current work, however,
we employ the same resolution of the two grids.
C
D
4
1
(a)
3
A
B
2
(b)
Figure 4: (a) Initial computation and material grids. (2) Deformed computation and
material grids. The red dots denote material particles and the black dots denote nodes
of the computation grid. The solid black lines show the connection between one material
particle (A) and computation nodes and the red dashed lines are the connections between
one computation node (3) and surrounding material particles.
This multi-grid strategy requires extra computation time. In the literature, one regular grid is used and the material grid is assumed to coincide
with the computation grid at all times [38, 40]. In the current work, we will
25
conduct simulations primarily using this single-grid simplification. Comparison with the multi-grid results are demonstrated. It will be shown that for
the examples considered, the mechanical responses computed from the multior single-grid approaches do not vary significantly although the deformed microstructures have distinct geometry.
The orientations of crystals also evolve with deformation. Upon convergence, the local crystallographic lattice rotations can be updated by the spin
rate calculated by
˜
ω̇(x) = Ω̇ + ω̇(x)
− ω̇ slip (x),
(63)
where Ω̇ = antisym(∇V) is the average spin rate over the microstructure
˜
˜ e (x) + ω̇
˜ p (x) is the spin rate fluctuation induced by elastic
domain. ω̇(x)
= ω̇
and plastic rotation. The last term that is subtracted is the rotation rate due
to plastic shear (slip) that does not distort the crystal lattice. It is calculated
by
ω̇ slip (x) =
Ns
X
β (α) · γ̇ (α) ,
(64)
α
where β (α) is the anti-symmetric Schmid tensor (β (α) = antisym(S(α) ) =
1 (α)
(s ⊗ n(α) − n(α) ⊗ s(α) ), and s(α) , n(α) are slip direction and normal to the
2
slip plane of the α-th slip system, respectively).
5. Crystal Plasticity Finite Element Method
We will compare the CEPFFT results with conventional full-field crystal
plasticity finite element method adopting a total Lagrangian scheme [27]
developed in our previous work [29]. The reference configuration of each
time step is the original undeformed microstructure. The total deformation
gradient F is multiplicatively decomposed into an elastic part Fe and a plastic
part Fp .
F = Fe Fp ,
(65)
where plastic deformation is volume conserved: |Fp | = 1.
The velocity gradient of plastic deformation Lp is determined by slip rates
on all slip systems:
X
Lp = Ḟp Fp−1 =
γ̇ (α) s(α) ⊗ n(α) .
(66)
α
26
The plastic strain rate and spin rate are calculated by ε̇p = sym(Lp ) and
ω̇ p = antisym(Lp ), respectively. The same flow rule (e.g. Eq. (39)) and
hardening law (e.g. Eq. (40)) as in the CEPFFT model are employed for
estimating shearing rates γ̇ (α) .
The orientation computed at Gauss points is updated using the elastic
deformation gradient as:
(α)
s(α) = Re s0 ,
(α)
n(α) = Re−T n0 ,
(α)
(67)
(α)
where Re is the rotation part of Fe , and s0 and n0 are slip direction and
plane normal of the α-th slip system of the undeformed grain.
6. Numerical examples
In this section, numerical examples conducted using the fast Fourier transform based approach are presented. Comparison between different formulations (e.g. visco-plasticity vs. elasto-plasticity) and between different methods (finite element method vs. fast Fourier transform based methods) are
conducted to validate the current developments. In addition, the basic FFT
methodology presented earlier is compared with the augmented Lagrangian
approach introduced in [40, 33]. Mechanical response and fatigue properties measured by strain based fatigue indicator parameters [50] of IN100
microstructures are studied using the novel CEPFFT method. The computational efficiency of CEPFFT and crystal plasticity finite element method
are presented to show the merit of the CEPFFT method. The use of the
multi-grid strategy is also discussed.
6.1. Basic formulation versus the augmented Lagrangian formulation
We will start with a benchmark plane strain example of 3D polycrystalline microstructure simulated using the crystal visco-plasticity fast Fourier
transform method (elastic response is neglected). The crystal visco-plastic
constitutive model is implemented in the basic framework highlighted earlier in this paper as well as using the so called augmented Lagrangian formulation [40, 33]. The single-grid strategy is adopted so that deformation
fluctuations are neglected.
FCC aluminum is considered. A cubic polycrystalline microstructure
composed of 64 grains is generated using the Voronoi tessellation scheme in
27
an 1mm3 domain. The microstructure is discretized by 16 × 16 × 16 equally
spaced voxels. The macroscale velocity gradient is


0.0 0.0 0.0
L = ∇V =  0.0 1.0 0.0  × 10−3 (s−1 ).
(68)
0.0 0.0 −1.0
The rate-dependent flow rule (Eq. (39)) is used with γ̇0 = 1s−1 and
m = 0.1. The hardening model given in Eq. (40) is adopted to update
slip resistances. The parameters in the hardening law are selected according
to [40]: κ0 = 47.0MP a, κ1 = 86.0MP a, θ0 = 550.0Mpa, and θ1 = 16.0MP a.
An arbitrary random texture is assigned to the microstructure. The initial
microstructure configuration and pole figures showing the random orientation
distribution are shown in Fig. 5.
Y
Y
1
0.8
X
X
X3
0.6
0.4
0.2
< 110 >
0
0
< 111 >
0
0.2
0.2
0.4
X2
0.4
0.6
0.6
0.8
X1
0.8
1
1
(a)
(b)
Figure 5: (a) The image representation of a 3D polycrystalline microstructure containing
64 grains. (b) Pole figures of the microstructure with randomly assigned orientations.
This example is used to validate the implementation of the general FFT
framework by providing a comparison with the alternative augmented Lagrangian implementation. After the thickness of the microstructure reduces
by 50%, the deformed microstructure and its stress and strain fields are plotted in Fig. 6. For both implementations, the stress distribution over the
microstructure is consistent with the grain geometry. The local mechanical responses of the two simulations are very close. Both intergranular and
intragranular heterogeneities of the stress and strain rate fields are captured.
The volume-averaged effective stress-strain responses of the entire microstructure computed by the two different algorithms are shown in Fig. 7(a).
28
Figure 6: Contour plots of plane strain deformed microstructures evaluated by different
algorithms. The top layer is the equivalent (plastic) strain field, and the bottom layer
is the equivalent stress field. (a) Crystal visco-plasticity fast Fourier transform approach
implemented in the basic formulation (b) Crystal visco-plasticity fast Fourier transform
approach implemented in the augmented Lagrangian formulation.
The two curves overlap. Pole figures showing the orientation distribution of
deformed microstructure are depicted in Fig. 7(b), from which we observe
the typical plane strain deformation texture pattern for both cases.
From the above comparison of the local and effective mechanical responses, we conclude that consistent results are obtained from the two algorithms. An error analysis study for this problem also has revealed that
the equilibrium error has similar convergence for both approaches. This fact
proves that there is no need for an augmented Lagrangian implementation of
the FFT algorithm and that the basic formulation has sufficiently accurate
performance for solving polycrystal plasticity problems [51].
It is also worth mentioning that the equilibrium error is mostly determined by the resolution of the microstructure. For high resolution, the equilibrium condition is fulfilled with smaller error (see Fig. 8). When the number
of pixels per side is doubled, the equilibrium error is approximately halved.
29
Y
250
X
Equivalent stress (MPa)
200
150
< 110 >
100
Y
X
< 111 >
Augmented
Lagrangian
Basic
Augmented
Lagrangian
50
0
0.1
0.2
0.3
0.4
0.5
Y
X
X
< 111 >
< 110 >
0
Y
Basic
0.6
Equivalent plastic strain
(a)
(b)
Figure 7: (a) The macroscopic effective stress-strain responses computed by the basic and
augmented Lagrangian crystal visco-plasticity FFT algorithms. (b) Pole figures of the
deformed microstructure texture.
6.2. Crystal elasto-plastic FFT simulations for polycrystalline microstructures
We consider here the same example as in Section 6.1 using the CEPFFT
method. Both the microstructure configuration and material parameters
are identical to the previous example. The elastic constants in the elastoplastic model are chosen to be C11 = 110 × 103 MP a, C12 = 59 × 103 MP a,
C44 = 26 × 103 MP a. The CEPFFT results are compared with the pure
visco-plastic computation as well as the results obtained from the crystal
plasticity finite element method. Both multi-grid and single-grid strategies
are also adopted and compared.
In the finite element simulation, homogeneous boundary condition is applied to the microstructure to drive its deformation while the boundary conditions of FFT-based simulations are periodic. The homogeneous boundary
condition enforces all boundary nodes to have the same deformation/velocity
gradient (e.g. Eq. (68) for the current plane strain problem), but heterogeneous nodal response inside the microstructure is allowed. This homogeneous
boundary condition will result in different local mechanical responses on the
microstructure from those obtained using the periodic boundary condition,
but the macroscopic effective response should be comparable. The FEM
simulation is conducted using our in-house solver extended based on [29].
We first adopt the main CEPFFT proposed in Section 2.4 with single-grid
strategy to perform simulations. The total strain, plastic strain, and stress
30
0.05
Pixel number Equilibrium
per side
error
Equilibrium error
0.04
0.03
8
4.292371e-02
16
2.195340e-02
32
1.116329e-02
64
5.609229e-03
0.02
0.01
0
0
10
20
30
40
50
60
70
Pixel number per side
Figure 8: Equilibrium error as a function of resolution (number of pixels per side). The
equilibrium error is evaluated when the convergence error reaches below 10−7 .
fields of the deformed microstructure after 50% thickness reduction computed
by different models are plotted in Fig. 9. It is seen that CEPFFT gives very
consistent prediction to both strain and stress fields with the pure viscoplastic results obtained in Section 6.1. The magnitude of the plastic strain
field prediction by CEPFFT is slightly smaller than the one predicted by
visco-plastic approach. On the other hand, the fields predicted by the FFTbased simulations show differences (especially for strain fields) from those
based on crystal plasticity finite element simulations. The major causes of the
differences between the FFT-based methods and the FEM approach include:
(1) the different boundary conditions (periodic by the FFT-based methods
and homogeneous by FEM), (2) the grain in FFT-based simulation is assigned
to the digital model in the unit of “point” while it is assigned to the finite
element model in the unit of “element” (each element contains 8 points),
(3) the fast Fourier transform results are computed directly on the voxel
points while the finite element results are computed on integration points
and extrapolated to nodes using a least squares method, and (4) algorithmic
differences. However, we notice that the stress fields exhibit similar patterns.
This is because the mean strains predicted by the two methods are about
the same and much larger in value than the strain fluctuations. One thing is
also worth mentioning is that the spatial variation of strain fields predicted
by crystal plasticity finite element method is actually on the same level with
CEPFFT. In the contour plots, the values shown in the finite element results
is adjusted to be smaller than in CEPFFT in order to show the stain variation
on boundary nodes, where homogeneous deformation is applied.
31
Figure 9: Contour plots of plane strain deformed microstructures evaluated by different
methods. The first row is the equivalent total strain field, the second row is the equivalent
plastic strain field, and the bottom row is the equivalent stress field. (a) Crystal viscoplasticity fast Fourier transform method. The total strain and plastic strain are identical
here since the elastic response is ignored in this model. (b) Crystal elasto-plasticity fast
Fourier transform method (CEPFFT) (c) Crystal plasticity finite element method.
The macroscopic effective stress-strain responses of the entire microstructure are compared in Fig. 10. The elastic response is successfully captured
by the CEPFFT model and is comparable to the finite element prediction. It
is observed that CEPFFT gives close prediction to the macroscopic response
with finite element approach even though the boundary conditions for the
two approaches are different.
The textures predicted by the three models are plotted using pole figures
in Fig. 11. It is observed that all the three sets of pole figures predicted are
very similar. The typical plane strain pattern is obtained.
To check the convergence of the results with respect to increasing resolution, we performed the same simulation on the same microstructure but
discretized by 32 × 32 × 32 voxels using CEPFFT. The local and effective
mechanical responses are compared in Figs. 12 and 13, respectively. It is observed that the two sets of results are consistent. Plotted in the two figures
are also the predictions using the modified CEPFFT approach introduced in
32
Figure 10: The macroscopic effective stress-strain responses of plane strain deformed microstructures predicted by different models. (a) Effective stress-total strain responses by
CEPFFT and crystal plasticity FEM; (b) The effective stress-plastic strain responses by
the three methods. Note that here CVPFFT denotes crystal visco-plasticity fast Fourier
method, and CPFEM refers to crystal plasticity finite element method.
Remark 2, where stress is linked to total strain rate through an elasto-plastic
modulus. We observe that the results from this modified implementation are
close to the main formulation of this paper (that computes elastic and plastic
fluctuations separately), although with milder spatial variation, implying its
good performance for this problem.
We also studied the equilibrium error of the main CEPFFT approach
for microstructures with different resolution. The improvement of error with
refining the image is seen from Fig. 14(a), which is similar to the viscoplasticity case. The convergence of the CEPFFT model as a function of
iteration number is shown in Fig. 14(b). We observe a fast convergence rate
for all tests.
In order to check the effect of heterogeneous deformation on mechanical
responses, we next repeat the above simulations with the multi-grid strategy. implemented with the CEPFFT. The deformed material grid obtained
by the multi-grid strategy and the homogeneous approximation are plotted
in Fig. 15. Contour plots of local mechanical responses are shown in Fig. 16.
Comparing with the single-grid results, we find that the predicted mechanical
responses do not vary much, although the deformed microstructures becomes
irregular. The effective stress-strain curve and grain orientation distribution
predicted by the multi-grid strategy are almost identical with the correspond33
Y
CVPFFT
Y
X
Y
X
< 111 >
< 110 >
(a)
CEPFFT
Y
X
Y
X
< 111 >
< 110 >
(b)
CPFEM
Y
X
X
< 111 >
< 110 >
(c)
Figure 11: Crystallographic textures, represented in pole figures, of plane strain deformed microstructures predicted by three models. (a) Crystal visco-plasticity fast Fourier
transform method (CVPFFT) (b) Crystal elasto-plasticity fast Fourier transform method
(CEPFFT) (c) Crystal plasticity finite element method (CPFEM).
ing results obtained with the single-grid prediction and are not repeated here.
We next examine an example of volume-conserved compression using
three different models. The same polycrystalline microstructure containing 64 grains and material parameters are used, while the imposed velocity
gradient becomes


0.5 0.0 0.0
L = ∇V =  0.0 0.5 0.0  × 10−3 (s−1 ).
(69)
0.0 0.0 −1.0
The single-grid strategy along with the main CEPFFT formulation is
firstly adopted. After the compression in z direction reaches 50%, the total
strain, plastic strain, and stress fields are plotted in Fig. 17. The local mechanical fields obtained from visco-plasticity and elasto-plasticity approaches
agree very well. We also observe that the finite element stress and strain field
results show very similar pattern with the FFT-based results despite the different kinematic conditions imposed on the FEM and FFT methodologies.
The macroscopic effective stress-strain responses of the entire microstructure are compared in Fig. 18. The elastic response is also successfully captured by the CEPFFT model and is comparable with the finite element
prediction.
The crystallographic textures described by pole figures are also captured
(Fig. 19). We can see once more that the different models give very close
prediction on the texture evolution of the volume-conserved compression.
To check the convergence of the results with respect to increasing resolution for volume-conserved compression deformation, we also performed the
34
Figure 12: Contour plots of plane strain deformed microstructures with different resolution
and methods. The first row is the equivalent total strain field, the second row is the
equivalent plastic strain field, and the bottom row is the equivalent stress field. (a) The
main CEPFFT method using 16 × 16 × 16-voxel microstructure. (b) The main CEPFFT
method using 32 × 32 × 32-voxel microstructure. (c) The modified CEPFFT (Remark 2)
using 16 × 16 × 16-voxel microstructure.
same simulation on the same microstructure but discretized by 32 × 32 × 32
voxels. The local and effective mechanical responses are compared in Figs. 20
and 21, respectively. The results computed using the modified CEPFFT formulation discussed in Remark 2 are also demonstrated. Both the local and
effective responses are close for the two implementations of the CEPFFT
method.
Comparing the multi-grid with the single-grid results, we note that the
irregular deformation is captured (Fig. 22), while the mechanical responses do
not exhibit obvious change (Fig. 23). Almost identical effective stress-strain
response and texture evolution with the single-grid prediction are observed
and as such these results are not shown here.
6.3. Investigation of fatigue indicator parameters of IN100
In this subsection, we adapt the CEPFFT method to study fatigue properties of Ni-based superalloys using the flow rule and hardening law for the
35
250
Y
Main CEPFFT-32PY
200
Equivalent stress (MPa)
X
X
80
< 111 >
< 110 >
Equivalent stress (MPa)
150
100
60
Modified CEPFFT-16P
40
Y
Main CEPFFT-16P
20
X
Modified CEPFFT-16P
50
0
Y
Main CEPFFT-32P
0
0.001
0.002
0.003
0.004
X
0.005
Equivalent strain
0
< 111 >
< 110 >
0
0.1
0.2
0.3
0.4
0.5
0.6
Equivalent strain
(a)
(b)
Figure 13: (a) The macroscopic effective stress-total strain responses by 16 × 16 × 16
microstructure and 32 × 32 × 32 microstructure obtained using different formulations.
(b) Crystallographic textures represented in pole figures. Main CEPFFT refers to the
main crystal elasto-plasticity FFT method and Modified CEPFFT refers to the crystal elasto-plasticity implementation using a homogeneous elasto-plastic medium approach
(Remark 2).
IN100 model introduced in Section 3.1.2. The same microstructure configuration as in the previous examples is used and the initial texture is random. The same parameters of the constitutive equations as listed in [30]
are used. The volume fractions and sizes of γ ′ precipitates are given by
fp1 = 0, fp2 = 0.42, d2 = 108nm, fp3 = 0.11, d3 = 7nm. The microstructure
is subjected to a 3-loop cyclic loading (tension and compression along the
z-direction). The stress-strain response in the z-direction during the 3 loading loops is plotted in Fig. 24 with comparison to an FEM simulation that
adopts the same constitutive model. The CEPFFT models implemented in
both the main formulation and the modified formulation (in Remark 2) are
tested. The single-grid update strategy is adopted. The microstructure input to the CEPFFT simulation is discretized by 16 × 16 × 16 voxels to be
consistent with the finite element input (16 ×16 ×16 cubic elements). Homogeneous boundary condition is applied to the finite element microstructure.
The CEPFFT simulations give similar prediction to the stress-strain “loop”
with the finite element model. Note that the loop predicted by the main
CEPFFT algorithm is wider than the loops obtained from the FEM or the
modified CEPFFT algorithm.
36
100
0.0014
Pixel number
per side
Equilibrium
error
10-1
8
1.339737e-03
10-2
16
6.971389e-04
32
3.562514e-04
64
1.789642e-04
8x8x8
16x16x16
32x32x32
64x64x64
0.001
0.0008
Convergence error
Equilibrium error
0.0012
10-3
10-4
0.0006
10-5
0.0004
10-6
0.0002
10-7
0
10
0
10
20
30
40
50
60
70
-8
1
1.5
2
2.5
Pixel number per side
Iterations
(a)
(b)
3
3.5
4
Figure 14: (a) Equilibrium error as a function of resolution (number of pixels per side).
The equilibrium error is evaluated using the basic formulation when the convergence error
reaches below 10−7 . (b) Convergence error versus number of iterations.
Figure 15: (a) Deformed microstructure predicted by multi-grid CEPFFT. (c) Deformed
microstructure predicted by single-grid CEPFFT.
We next compute strain based fatigue indicator parameters (FIPs) related
to small crack formation and early growth [50]. The three fatigue indicator parameters of interest are the cumulative plastic strain per cycle (Pcyc ),
which correlates to the crack incubation life; the Fatemi-Socie parameter
(PF S ), which relates to the small crack growth; and the maximum range of
cyclic plastic shear strain parameter (Pmps ). The definitions of these fatigue
indicator parameters can be found in [44]. Contour plots of fatigue indicator
parameters over the microstructure are plotted in Fig. 25. All fatigue indicator parameters are computed for the last (3rd) loading loop. We observe that
the fatigue indicator parameters contour plots demonstrate similar patterns
in CEPFFT and crystal plasticity finite element simulations. The results
37
Figure 16: Contour plots of plane strain deformed microstructures computed using different strategies. The first row shows results obtained from the multi-grid CEPFFT, and
the bottom row shows results obtained from the single-grid CEPFFT. (a) Equivalent total
strain (b) Equivalent plastic strain (c) Equivalent stress.
predicted by the main CEPFFT approach show more heterogeneity and are
closer to finite element results than the modified FFT model (Remark 2).
The convergence of the CEPFFT simulation with respect to resolution
and the multi-grid effect are also tested for this fatigue problem. The fatigue indicator parameters are plotted in Fig. 26. It is observed that the
deformation heterogeneity does not affect significantly the fatigue indicator
parameters fields and the distortion of the microstructure is small since the
total strain is small for this problem. The finer microstructure (32 × 32 × 32
voxels) gives consistent prediction to local fatigue indicator parameters with
the coarse grid prediction, while larger heterogeneity is captured.
Grain level fatigue indicator parameters, namely the maximum and average fatigue indicator parameters of voxels within individual grains [44] are
also extracted. Their distributions are shown in Fig. 27. The distributions
computed using the main CEPFFT based on the 16 × 16 × 16-voxel microstructure are shown in (a). The ones computed using the main CEPFFT
on the 32 ×32 ×32-voxel microstructure are plotted in (b). The results of the
modified implementation of CEPFFT (Remark 2) are depicted in (c), and
the crystal plasticity finite element results are given in (d). These results are
obtained using the single-grid strategy. In the figure, the FIP with “max”
prefix means the maximum value of the FIP of voxels within an individual
38
Figure 17: Contour plots of compressed microstructures evaluated by different methods.
The first row is equivalent total strain field, the second row is equivalent plastic strain
field, and the bottom row is equivalent stress field. (a) Crystal visco-plasticity fast Fourier
transform method (CVPFFT). The total strain and plastic strain are identical here since
the elastic response is ignored in this model. (b) Crystal elasto-plasticity fast Fourier
transform method (CEPFFT). (c) Crystal plasticity finite element method (CPFEM).
grain. The one with “ave” prefix indicates that it is the average FIP over
voxels within an individual grain. The distributions are normalized by the
maximum value of each FIP over all grains of the microstructure. It is seen
that the distributions of grain level fatigue indicator parameters predicted by
the main CEPFFT on the finer microstructure are closer to the finite element
results. The ones predicted by the main CEPFFT on coarse microstructure
are also close to the fine microstructure estimations. On the other hand, the
maximum and average fatigue indicator parameters distributions predicted
by the modified CEPFFT are not very distinguishable, implying that the
intragranular heterogeneity of fatigue indicator parameters predicted by the
modified CEPFFT formulation is weak. This is consistent with the contour
plots of the FIP fields. A more sophisticated elasto-plastic modulus may
resolve this problem.
39
Figure 18: The macroscopic effective stress-strain responses of compressed microstructures
predicted by different models. (a) Effective stress-total strain responses by CEPFFT
and crystal plasticity finite element method (CPFEM); (b) Effective stress-plastic strain
responses by all three methods. CVPFFT refers to crystal visco-plasticity fast Fourier
method.
6.4. Computation Efficiency
The computational efficiency of CEPFFT method is important since our
goal is to develop a high efficient full-field simulator that can be adopted in
stochastic simulations as the deterministic solver to further study the probabilistic nature of material properties. As mentioned before, with the employment of fast Fourier transform to solve for local mechanical responses, the
cumbersome matrix inversion in finite element simulation is circumvented,
which leads to great improvement of the computation speed. The computation times of microstructures subjected to plane strain deformation to 0.02
strain are shown in Fig. 28(a) and tabulated in Table 1. All tests were run on
the Teragrid TACC Lonestar Linux Cluster. Each computation node contains
two Xeon Intel Hexa-Core 64-bit Westmere processors, the core frequency of
which is 3.33GHz. Single-grid morphology evolution strategy and main formulation are adopted by all CEPFFT tests. We first ran examples on one
single processor, respectively. The CEPFFT simulation of a microstructure
discretized by 163 voxels only takes less than 3 minutes, while the finite element simulation of a microstructure containing 163 elements takes 13 hours.
With an increased resolution of 323 voxels, the computation time of CEPFFT
increases to about 22 minutes, which is still significantly smaller than the finite element simulation. In order to further reduce the computation time,
40
Y
CVPFFT
Y
X
Y
X
< 111 >
< 110 >
(a)
CEPFFT
Y
X
Y
X
< 111 >
< 110 >
(b)
CPFEM
Y
X
X
< 111 >
< 110 >
(c)
Figure 19: Crystallographic textures, represented in pole figures, of compressed microstructures predicted by three models. (a) Crystal visco-plasticity fast Fourier transform method (CVPFFT). (b) Crystal elasto-plasticity fast Fourier transform method
(CEPFFT). (c) Crystal plasticity finite element method (CPFEM).
Table 1: Computation times for microstructures under plane strain and cyclic deformations simulated using different methods. CPFEM refers to crystal plasticity finite element
method.
1-processor, 60-processor, 1-processor, 240-processor, 1-processor, 240-processor,
CEPFFT, CEPFFT, CEPFFT,
CEPFFT,
CPFEM,
CPFEM,
163 -voxel
163 -voxel
323 -voxel
323 -voxel 163 -element 163 -element
Plane strain
164s
5s
1309s
17s
46389s
560s
Cyclic loading
554s
16s
4385s
57s
46873s
522s
we also parallelized our code. The 163 -voxel CEPFFT simulation only takes
5 seconds to finish if run on 60 processors. The 323 -voxel simulation running
on 240 processors takes 17 seconds to finish. On the other hand, the finite
element simulation of 163 -element microstructure takes about 10 minutes to
complete when 240 processors are used.
Besides plane strain deformation, we also tested the computational efficiency of the microstructure fatigue problem. The computation times of
one complete loading loop are captured and plotted in Fig. 28(b). We observe, again, that the CEPFFT simulation is much faster than the crystal
plasticity finite element method. Note that the time increment required by
the finite element simulation is smaller than that required by the fast Fourier
transform simulation in order to reach convergence. However, even if the two
approaches use the same time increment, the FFT computation is still much
faster than the FEM simulation for the microstructure resolution.
With the adoption of fast Fourier transform for solving the polycrystalline
microstructure boundary value problem, the most time consuming operation
41
Figure 20: Contour plots of compressed microstructures with different resolution and
formulation. The first row is the equivalent total strain field, the second row is the equivalent plastic strain field, and the bottom row is the equivalent stress field. (a) The main
CEPFFT using 16 × 16 × 16-voxel microstructure. (b) The main CEPFFT method using
32 × 32 × 32-voxel microstructure (c) The CEPFFT method implemented using a homogeneous elasto-plastic medium approach (Remark 2) using 16 × 16 × 16-voxel microstructure.
in the FFT-based simulation is then the nonlinear constitutive calculations
defined by the crystal plasticity theory. Another branch of research using
Fourier methods in computational crystal plasticity has been conducted by
Kalidindi and coworkers, who have introduced Fourier expansion to establish
hypersurfaces of constitutive equations so that the iterative computation in
updating nonlinear local properties can be approximated using spectral interpolation [52, 53, 54, 55]. Different from our strategy of solving the overall microstructure governing equations using Fourier transform, they replace
nonlinear local constitutive equations by Fourier expansion after calibrating
coefficients of the Fourier series using a given database. The computation
efficiency of a crystal plasticity problem may be further improved by coupling their spectral method with the present crystal plasticity fast Fourier
transform framework.
42
250
Y
Main CEPFFT-32P
200
Y
Equivalent stress (MPa)
X
X
80
Equivalent stress (MPa)
100
60
Modified CEPFFT-16P
40
Y
Main CEPFFT-16P
X
Modified CEPFFT-16P
0
0
0.001
0.002
0.003
0.004
0.1
0.2
0.3
X
0.005
Equivalent strain
0
Y
Main CEPFFT-32P
20
50
0
< 111 >
< 110 >
150
< 111 >
< 110 >
0.4
0.5
Equivalent strain
(a)
(b)
Figure 21: (a) The macroscopic effective stress-total strain responses by 16 × 16 × 16
microstructure and 32 × 32 × 32 microstructure obtained using different formulations.
(b) Crystallographic textures represented in pole figures. Main CEPFFT refers to the
main crystal elasto-plasticity FFT method and Modified CEPFFT refers to the crystal elasto-plasticity implementation using a homogeneous elasto-plastic medium approach
(Remark 2).
7. Conclusions
In this work, we developed an efficient FFT full field model to investigate elasto-plastic properties of polycrystalline materials by interrogating
image-based realistic microstructures. The elastic and plastic responses were
computed separately. An integrated formulation was also proposed using a
particular choice for the elasto-plastic moduli. Predictive capability and computational efficiency of the newly developed CEPFFT method were presented
using numerical examples in comparison with visco-plasticity approach and
crystal plasticity finite element simulations. Error tests were conducted to
show the comparable performance of the basic and augmented Lagrangian
algorithms. The multi-grid strategy was implemented to allow comparison
with the single-grid simplification. The self-consistence (convergence with
refining discretization) of the CEPFFT method was shown through simulations with increasing resolution of the microstructure. Fatigue properties of
Ni-based superalloy microstructures described by fatigue indicator parameters were studied using the CEPFFT method. The main discoveries of this
work are as follows:
43
Figure 22: (a) Deformed microstructure predicted by multi-grid CEPFFT. (c) Deformed
microstructure predicted by single-grid CEPFFT.
1. The elastic response of the microstructure is successfully captured by
CEPFFT. Fatigue properties of nickel-based superalloy microstructures
can be efficiently investigated using CEPFFT method. The computed
results are comparable with the those based on crystal plasticity finite
element simulation.
2. The equilibrium error resulted from the augmented Lagrangian algorithm is comparable with the one resulting from the basic formulation
and depends on the resolution of the input image.
3. The CEPFFT model provides consistent prediction of the macroscopic
effective mechanical responses and local stress field with the finite element and crystal visco-plasticity fast Fourier transform simulations.
The crystallographic texture patterns of microstructures under different
deformation modes estimated by the three methods are almost identical.
4. The local strain fields obtained by the CEPFFT implementation agree
very well with visco-plastic results and have similar patterns with the
crystal plasticity finite element results.
5. The multi-grid strategy that uses separate computation and material grids was shown to successfully capture irregular configuration of
deformed microstructures. However, it does not result in significant
changes in the mechanical response for the tested examples in comparison to the single-grid method.
6. The required computation time of CEPFFT is significantly less than
that of crystal plasticity finite element method since the methodology
does not require the solution of any matrix equations as in the FEM.
44
Figure 23: Contour plots of compressed microstructures computed using the multi-grid
strategy. The first row shows results obtained from the multi-grid CEPFFT, and the
bottom row shows results obtained from the single-grid CEPFFT. (a) Equivalent total
strain. (b) Equivalent plastic strain. (c) Equivalent stress.
This computational efficiency provides significant potential in integrating the FFT approach with stochastic multiscale materials simulations.
7. Various CEPFFT can be introduced. The approach highlighted in this
paper computes separately the elastic and plastic responses. Careful design of the elasto-plastic modulus for the homogeneous reference
medium can improve the computational efficiency.
Acknowledgements
This research was supported by an OSD/AFOSR MURI09 award on uncertainty quantification, the U.S. Department of Energy, Office of Science,
Advanced Scientific Computing Research and the Computational Mathematics program of the National Science Foundation (NSF) (award DMS0809062). This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the
U.S. Department of Energy under Contract No. DE-AC02-05CH11231. Additional computing resources were provided by the NSF through TeraGrid
resources provided by NCSA under grant number TG-DMS090007.
45
Stress
500
0
Main CEPFFT
Modified CEPFFT
CPFEM
-500
-0.006
-0.004
-0.002
0
0.002
0.004
0.006
Strain
Figure 24: Stress-strain responses of three loops of cyclic loading computed by CEPFFT
and CPFEM. Main CEPFFT refers to the main crystal elasto-plasticity FFT method,
modified CEPFFT refers to the crystal elasto-plasticity implementation using a homogeneous elasto-plastic medium approach (Remark 2), and CPFEM is the crystal plasticity
finite element method.
Figure 25: Contour plots of fatigue indicator parameters fields. Main CEPFFT results
are placed in the top row, results from the modified CEPFFT formulation (Remark 2) are
placed in the middle row, and crystal plasticity finite element results are located in the
bottom row. (a) Pcyc , (b) PF S , (c) Pmps .
46
Figure 26: Contour plots of fatigue indicator parameters fields evaluated by the multi-grid
CEPFFT on coarse microstructure (top row), single-grid CEPFFT on coarse microstructure (middle row), and single-grid CEPFFT on fine microstructure (bottom row). (a)
Pcyc , (b) PF S , (c) Pmps .
47
0.3
0.25
0.2
maxPcyc
maxPFS
0.15
maxPmps
avePcyc
0.1
avePFS
0.05
Fraction of grains
Fraction of grains
0.25
avePmps
maxPcyc
0.2
maxPFS
maxPmps
0.15
avePcyc
0.1
avePFS
avePmps
0.05
0
0
0
0.2
0.4
0.6
0.8
Grain level FIP/max FIP
0
1
0.2
0.4
0.6
0.8
1
Grain level FIP/max FIP
(a)
(b)
0.35
0.25
maxPcyc
maxPFS
0.15
maxPmps
avePcyc
0.1
avePFS
0.05
Fraction of grains
ains
Fraction of grains
ains
0.3
0.2
avePmps
maxPcyc
0.25
maxPFS
0.2
maxPmps
0.15
avePcyc
0.1
avePFS
avePmps
0.05
0
0
0
0.2
0.4
0.6
0.8
0
1
0.2
Grain level FIP/max FIP
0.4
0.6
0.8
1
Grain level FIP/max FIP
(c)
(d)
Figure 27: Distribution of the fatigue indicator parameters among grains computed by (a)
Main CEPFFT with coarse microstructure, (b) Main CEPFFT with fine microstructure,
(c) Modified CEPFFT (Remark 2) with coarse microstructure and (d) crystal plasticity
finite element method. These distributions are normalized by the maximum values of each
FIP over all grains.
100000
Computation time (s)
Computation time (s)
100000
10000
1000
100
10
1
(a)
10000
1000
100
10
1
(b)
Figure 28: Computation times of simulations using different methods on microstructures
with different resolutions. (a) Plane strain deformation to 0.02 strain. (b) One complete
loop of cyclic loading on the IN100 microstructure. The computation times are shown in
logarithmic scale. In the figure, 1ProcCEPFFT-16P means the CEPFFT simulation of a
microstructure discretized by 16 × 16 × 16 voxels using 1 processor and 1ProcCPFEM-16E
means the crystal plasticity finite element simulation of a microstructure discretized by
16 × 16 × 16 elements using 1 processor.
48
References
[1] G. Sachs, Plasticity problems in metals, Trans. Faraday Soc. 24 (1928)
84–92.
[2] G. Taylor, Plastic strain in metals, Journal of the Inistitute of Metals
62 (1938) 307–324.
[3] S. Ahzi, R. J. Asaro, D. M. Parks, Application of crystal plasticity
theory for mechanically processed bscco superconductors, Mechanics of
Materials 15 (1993) 201–222.
[4] T. Lin, Analysis of elastic and plastic strains of a face-centred cubic
crystal, Journal of the Mechanics and Physics of Solids 5 (1957) 143–
149.
[5] J. W. Hutchinson, Bounds and self-consistent estimates for creep of
polycrystalline materials, Proceedings of the Royal Society of London.
A. Mathematical and Physical Sciences 348 (1976) 101–127.
[6] R. Asaro, A. Needleman, Overview no. 42 texture development and
strain hardening in rate dependent polycrystals, Acta Metallurgica 33
(1985) 923–953.
[7] S. Kalidindi, C. Bronkhorst, L. Anand, Crystallographic texture evolution in bulk deformation processing of fcc metals, Journal of the Mechanics and Physics of Solids 40 (1992) 537–569.
[8] S. R. Kalidindi, Incorporation of deformation twinning in crystal plasticity models, Journal of the Mechanics and Physics of Solids 46 (1998)
267–290.
[9] P. R. Dawson, S. R. MacEwen, P.-D. Wu, Advances in sheet metal forming analyses: dealing with mechanical anisotropy from crystallographic
texture, International Materials Reviews 48 (2003) 86–122.
[10] B. Kouchmeshky, N. Zabaras, Modeling the response of hcp polycrystals
deforming by slip and twinning using a finite element representation of
the orientation space, Computational Materials Science 45 (2009) 1043–
1051.
49
[11] E. Kröner, On the plastic deformation of polycrystals, Acta Metallurgica
9 (1961) 155–161.
[12] J. D. Eshelby, The determination of the elastic field of an ellipsoidal
inclusion, and related problems, Proceedings of the Royal Society of
London. Series A, Mathematical and Physical Sciences 241 (1957) 376–
396.
[13] R. Hill, Continuum micro-mechanics of elastoplastic polycrystals, Journal of the Mechanics and Physics of Solids 13 (1965) 89–101.
[14] J. W. Hutchinson, Elastic-plastic behaviour of polycrystalline metals
and composites, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 319 (1970) 247–272.
[15] T. Iwakuma, S. Nemat-Nasser, Finite elastic-plastic deformation of
polycrystalline metals, Proceedings of the Royal Society of London.
A. Mathematical and Physical Sciences 394 (1984) 87–119.
[16] A. Molinari, G. Canova, S. Ahzi, A self consistent approach of the large
deformation polycrystal viscoplasticity, Acta Metallurgica 35 (1987)
2983–2994.
[17] R. Lebensohn, C. Tomé, A self-consistent anisotropic approach for the
simulation of plastic deformation and texture development of polycrystals: Application to zirconium alloys, Acta Metallurgica et Materialia
41 (1993) 2611–2624.
[18] R. Lebensohn, C. Tomé, A self-consistent viscoplastic model: prediction
of rolling textures of anisotropic polycrystals, Materials Science and
Engineering: A 175 (1994) 71–82. NATO Advanced Research Workshop
on Polyphase Polycrystal Plasticity.
[19] H. Wang, P. Wu, C. Tomé, Y. Huang, A finite strain elastic–viscoplastic
self-consistent model for polycrystalline materials, Journal of the Mechanics and Physics of Solids 58 (2010) 594–612.
[20] R. Masson, A. Zaoui, Self-consistent estimates for the rate-dependent
elastoplastic behaviour of polycrystalline materials, Journal of the Mechanics and Physics of Solids 47 (1999) 1543–1568.
50
[21] A. Zaoui, R. Masson, Micromechanics-based modeling of plastic polycrystals: an affine formulation, Materials Science and Engineering: A
285 (2000) 418–424.
[22] P. P. Castañeda, Second-order homogenization estimates for nonlinear
composites incorporating field fluctuations: I-theory, Journal of the
Mechanics and Physics of Solids 50 (2002) 737–757.
[23] Y. Liu, P. P. Castañeda, Second-order theory for the effective behavior and field fluctuations in viscoplastic polycrystals, Journal of the
Mechanics and Physics of Solids 52 (2004) 467–495.
[24] R. A. Lebensohn, C. N. Tomé, P. P. Castañeda, Self-consistent modelling
of the mechanical behaviour of viscoplastic polycrystals incorporating
intragranular field fluctuations, Philosophical Magazine 87 (2007) 4287–
4322.
[25] D. Peirce, R. Asaro, A. Needleman, An analysis of nonuniform and
localized deformation in ductile single crystals, Acta Metallurgica 30
(1982) 1087–1119.
[26] F. Roters, P. Eisenlohr, L. Hantcherli, D. Tjahjanto, T. Bieler, D. Raabe,
Overview of constitutive laws, kinematics, homogenization and multiscale methods in crystal plasticity finite-element modeling: Theory, experiments, applications, Acta Materialia 58 (2010) 1152–1211.
[27] L. Anand, M. Kothari, A computational procedure for rate-independent
crystal plasticity, Journal of the Mechanics and Physics of Solids 44
(1996) 525–558.
[28] A. Acharya, A. Beaudoin, Grain-size effect in viscoplastic polycrystals
at moderate strains, Journal of the Mechanics and Physics of Solids 48
(2000) 2213–2230.
[29] W. Li, N. Zabaras, A virtual environment for the interrogation of 3d
polycrystalline microstructures including grain size effects, Computational Materials Science 44 (2009) 1163–1177.
[30] C. P. Przybyla, D. L. McDowell, Microstructure-sensitive extreme value
probabilities for high cycle fatigue of Ni-base superalloy IN100, International Journal of Plasticity 26 (2010) 372–394.
51
[31] B. Wen, N. Zabaras, Investigating variability of fatigue indicator parameters of two-phase nickel-based superalloy microstructures, Computational Materials Science 51 (2012) 455–481.
[32] R. Lebensohn, N-site modeling of a 3D viscoplastic polycrystal using
Fast Fourier Transform, Acta Materialia 49 (2001) 2723–2737.
[33] R. Lebensohn, P. Ponte Castañeda, R. Brenner, O. Castelnau, Fullfield vs. homogenization methods to predict microstructure-property
relationships of polycrystalline materials, in: S. Ghosh, D. Dimiduk
(Eds.), Computational Methods for Microstructure-Property Relationships, Springer, 2011, pp. 393–441.
[34] H. Moulinec, P. Suquet, A fast numerical method for computing the
linear and nonlinear mechanical properties of composites, Comptes rendus de l’Académie des sciences. Série II, Mécanique, Physique, Chimie,
Astronomie 318 (1994) 1417–1423.
[35] H. Moulinec, P. Suquet, A numerical method for computing the overall
response of nonlinear composites with complex microstructure, Computer Methods in Applied Mechanics and Engineering 157 (1998) 69–94.
[36] J. C. Michel, H. Moulinec, P. Suquet, A Computational Method Based
on Augmented Lagrangians and Fast Fourier Transforms for Composites with High Contrast, CMES-Computer Modeling in Engineering &
Sciences 1 (2000) 79–88.
[37] R. Lebensohn, O. Castelnau, R. Brenner, P. Gilormini, Study of the antiplane deformation of linear 2-d polycrystals with different microstructures, International Journal of Solids and Structures 42 (2005) 5441–
5459.
[38] R. A. Lebensohn, R. Brenner, O. Castelnau, A. D. Rollett, Orientation
image-based micromechanical modelling of subgrain texture evolution
in polycrystalline copper, Acta Materialia 56 (2008) 3914–3926.
[39] R. A. Lebensohn, Y. Liu, P. P. Castañeda, Macroscopic properties and
field fluctuations in model power-law polycrystals: full-field solutions
versus self-consistent estimates, Proceedings of the Royal Society of
London. Series A: Mathematical, Physical and Engineering Sciences 460
(2004) 1381–1405.
52
[40] A. Prakash, R. A. Lebensohn, Simulation of micromechanical behavior
of polycrystals: finite elements versus fast fourier transforms, Modelling
and Simulation in Materials Science and Engineering 17 (2009) 064010.
[41] B. Liu, D. Raabe, F. Roters, P. Eisenlohr, R. A. Lebensohn, Comparison
of finite element and fast fourier transform crystal plasticity solvers for
texture prediction, Modelling and Simulation in Materials Science and
Engineering 18 (2010) 085005.
[42] S.-J. Kim, D. H. Kim, K. H. Oh, A. D. Rollett, R. A. Lebensohn, H. N.
Han, An elastoplastic finite element modeling coupled with orientation
image based micromechanical approach, AIP Conference Proceedings
1252 (2010) 103–106.
[43] N. Lahellec, J. Michel, H. Moulinec, P. Suquet, Analysis of inhomogeneous materials at large strains using fast Fourier transforms, in:
C. Miehe (Ed.), IUTAM symposium on computational mechanics of
solids materials, Kluwer Academic, Stuttgart, Germany, 2001, pp. 247–
268.
[44] M. Shenoy, J. Zhang, D. McDowell, Estimating fatigue sensitivity to
polycrystalline ni-base superalloy microstructures using a computational
approach, Fatigue & Fracture of Engineering Materials & Structures 30
(2007) 889–904.
[45] M. Shenoy, Y. Tjiptowidjojo, D. McDowell, Microstructure-sensitive
modeling of polycrystalline IN 100, International Journal of Plasticity
24 (2008) 1694–1730. Special Issue in Honor of Jean-Louis Chaboche.
[46] S. Sankaran, N. Zabaras, Computing property variability of polycrystals
induced by grain size and orientation uncertainties, Acta Materialia 55
(2007) 2279–2290.
[47] M. Frigo, S. Johnson, The design and implementation of FFTW3, Proceedings of the IEEE 93 (2005) 216–231.
[48] D. Sulsky, Z. Chen, H. Schreyer, A particle method for historydependent materials, Computer Methods in Applied Mechanics and
Engineering 118 (1994) 179 – 196.
53
[49] D. Sulsky, S.-J. Zhou, H. L. Schreyer, Application of a particle-incell method to solid mechanics, Computer Physics Communications 87
(1995) 236 – 252. Particle Simulation Methods.
[50] D. McDowell, F. Dunne, Microstructure-sensitive computational modeling of fatigue crack formation, International Journal of Fatigue 32
(2010) 1521–1542.
[51] P. Suquet, personal communication, 2011.
[52] S. R. Kalidindi, H. K. Duvvuru, Spectral methods for capturing crystallographic texture evolution during large plastic strains in metals, Acta
Materialia 53 (2005) 3613–3623.
[53] S. R. Kalidindi, H. K. Duvvuru, M. Knezevic, Spectral calibration of
crystal plasticity models, Acta Materialia 54 (2006) 1795–1804.
[54] M. Knezevic, S. R. Kalidindi, D. Fullwood, Computationally efficient
database and spectral interpolation for fully plastic taylor-type crystal
plasticity calculations of face-centered cubic polycrystals, International
Journal of Plasticity 24 (2008) 1264–1276.
[55] M. Knezevic, H. F. Al-Harbi, S. R. Kalidindi, Crystal plasticity simulations using discrete fourier transforms, Acta Materialia 57 (2009)
1777–1784.
54
Fly UP