PaperlayerView exact revision on arXiv

Source-linked AI summary

Finite strain FFT-based non-linear solvers made simple

T. W. J. de Geus, J. Vondrejc, J. Zeman, R. H. J. Peerlings, M. G. D. Geers

arXiv:1603.08893v3math.NAcs.NA

TL;DR

Periodic-cell homogenization needs efficient equilibrium solutions for complex microstructures, but nonlinear FFT methods have relied on a reference medium. This paper extends a variational FFT formulation to finite strains, using a projection operator, Newton iterations, and an iterative linear solver. Its examples demonstrate use with different material models, while the simulations illustrate local phase responses and convergence behavior.

  • Problem

    Nonlinear FFT-based finite-strain solvers retain an auxiliary reference medium, which may affect convergence and obscures consistent linearization.

    Source: §1 · p. 2

  • Method

    The paper extends a variational FFT formulation to finite strains, using a compatibility projection with Newton iterations and an iterative linear solver.

    Source: §1 · p. 2, §2.2 · p. 3

  • Results

    The two examples demonstrate the approach with different material models; the elasto-plastic example averages 2.2 iterations per increment over 250 increments.

    Source: §1 · p. 2, §5.1 · p. 10, §5.5 · p. 12

  • Takeaways & Limitations

    The examples show the formulation can treat finite-strain models defined in either the undeformed or deformed configuration, including an elasto-plastic model in the current configuration.

    Source: §1 · p. 2, §5.1 · p. 10

    The study restricts its boundary conditions to a fully prescribed macroscopic deformation gradient, identified as the easiest and most efficient choice.

    Source: §3.1 · p. 7

Abstract

from arXiv · show

Computational micromechanics and homogenization require the solution of the mechanical equilibrium of a periodic cell that comprises a (generally complex) microstructure. Techniques that apply the Fast Fourier Transform have attracted much attention as they outperform other methods in terms of speed and memory footprint. Moreover, the Fast Fourier Transform is a natural companion of pixel-based digital images which often serve as input. In its original form, one of the biggest challenges for the method is the treatment of (geometrically) non-linear problems, partially due to the need for a uniform linear reference problem. In a geometrically linear setting, the problem has recently been treated in a variational form resulting in an unconditionally stable scheme that combines Newton iterations with an iterative linear solver, and therefore exhibits robust and quadratic convergence behavior. Through this approach, well-known key ingredients were recovered in terms of discretization, numerical quadrature, consistent linearization of the material model, and the iterative solution of the resulting linear system. As a result, the extension to finite strains, using arbitrary constitutive models, is at hand. Because of the application of the Fast Fourier Transform, the implementation is substantially easier than that of other (Finite Element) methods. Both claims are demonstrated in this paper and substantiated with a simple code in Python of just 59 lines (without comments). The aim is to render the method transparent and accessible, whereby researchers that are new to this method should be able to implement it efficiently. The potential of this method is demonstrated using two examples, each with a different material model.

1 Introduction

The paper extends a variational FFT framework to finite-strain equilibrium, addressing the reference-medium dependence of earlier nonlinear formulations. It presents the method as straightforward to implement and demonstrates it with different constitutive models.

1 Introduction

  • FFT methods offer speed and memory-footprint advantages over finite elements for periodic cells, while pairing naturally with pixel-based microstructure representations.The paper identifies high-resolution, efficient three-dimensional representation of complex microstructures as a motivation for computational homogenization.Source: §1 · p. 1
  • Earlier nonlinear FFT extensions retained a reference medium, which may affect convergence and obscures consistent linearization.The paper notes that an available open-source implementation offers different solution techniques but does not include consistent linearization.Source: §1 · p. 2
  • The variational formulation removes the auxiliary elastic problem and separates discretization, quadrature, constitutive linearization, and linear-system solution.Its projection maps arbitrary tensor fields to compatible ones; coincident nodes and integration points support local strain and stress updates, while the local consistent tangent enables quadratic convergence.Source: §1 · p. 2
  • The paper extends the variational approach to finite strains and presents it as applicable to arbitrary constitutive models defined in either the undeformed or deformed configuration.For models defined in the deformed configuration, stress and consistent tangent are pulled back to the undeformed configuration using tensor operations.Source: §1 · p. 2
  • The paper includes a compact Python implementation and makes its code freely available to support independent implementation and further examples.The introduction describes the code as well below 200 lines and invites researchers to contribute additional examples.Source: §1 · p. 2
  • The paper extends prior variational finite-strain work, presents a pseudo-algorithm, and demonstrates the approach with two examples using different material models.It also points readers to the Python code and an appendix explaining the projection operator.Source: §1 · p. 2

2 Method

The method formulates periodic finite-strain equilibrium variationally, discretizes deformation gradients on a regular grid, and enforces compatibility with a Fourier-space projection. Newton iterations and an iterative linear solver then determine equilibrium without a constitutive-specific reference medium.

2.1 Problem statement

  • The problem is periodic-cell equilibrium for a prescribed overall deformation, with a regular pixel or voxel grid representing one or more material phases.The derivation uses nonlinear elasticity for notation, while the paper says extension to time- or history-dependent materials is straightforward and includes those in its examples.Source: §2.1 · p. 2

2.2 Weak form

  • The weak form expresses equilibrium using virtual deformation gradients and stress, with compatibility imposed as a constraint rather than built into the test functions.Unlike the finite-element position or displacement interpolation, the FFT formulation works in terms of strains or, for finite deformation, the deformation gradient.Source: §2.2 · p. 3

2.3 Projection to a compatible solution space

  • The projection operator G maps arbitrary fields to compatible fields, allowing the weak form to be tested with arbitrary periodic fields while the deformation gradient remains compatible.Its Fourier-space expression is closed-form, and convolution there becomes a local double tensor contraction.Source: §2.3 · p. 4

2.4 Discretization

  • A Galerkin scheme approximates the deformation-gradient field and test fields with the same nodal values and trigonometric shape functions.The shape functions are interpolatory and form a partition of unity.Source: §2.4 · p. 4

2.5 Quadrature

  • Trapezoidal quadrature gives equal weight to nodes and integration points, and the Fourier-space projection enforces a zero-mean field at zero frequency.The scaled frequency vector accounts for periodic-cell size; even-sized grids require additional treatment to enforce zero mean along with compatibility.Source: §2.5 · p. 5

2.6 Linearization

  • Newton iterations linearize the nodal equilibrium equations using local constitutive tangents, and an iterative solver enforces compatibility in the linear updates.The paper lists conjugate gradient, GMRES, Chebyshev, and Richardson iterations as options; an alternative method satisfies compatibility only at convergence.Source: §2.6 · p. 6

2.7 Connection to existing FFT-based solvers

  • Unlike reference-medium Green operators, G enforces compatibility through a constitutive-law-independent projection and remains constant throughout the simulation.Earlier approaches use an auxiliary constant stiffness, selected heuristically and requiring the Green operator to be updated at every Newton iteration, which can be costly.Source: §2.7 · p. 6

3 Implementation

The implementation solves equilibrium incrementally with Newton iterations and conjugate gradient linear solves. It prescribes the macroscopic deformation, provides a Python implementation, and avoids assembling the full system matrix.

3 Implementation

  • Each increment updates nodal deformation gradients through Newton iterations, with conjugate gradient solving the linearized system until the specified tolerances are met.The Newton and conjugate-gradient solves use separate accuracy criteria.Source: §3 · p. 7

3.1 Boundary conditions

  • Prescribing the macroscopic deformation gradient makes each iteration satisfy that deformation exactly, while the tangent distributes its increment across the heterogeneous microstructure.The paper identifies this treatment with essential (Dirichlet) boundary conditions in FEM.Source: §3.1 · p. 7

3.2 Program structure

  • The accompanying Python implementation uses standard scientific libraries, with no custom software or libraries required.The code is provided for the hyper-elastic example in Appendix B.Source: §3.2 · p. 7

Algorithm 1

  • Conjugate gradient repeatedly applies the operator to a deformation increment; storing its components separately avoids assembling the fully populated system matrix.Source: Algorithm 1 · p. 9
  • In three dimensions, second-order tensor fields are stored as five-dimensional arrays and fourth-order fields as seven-dimensional arrays, without requiring sparse-matrix representations.The passage gives storage sizes of 3^2 × n and 3^4 × n entries, respectively.Source: Algorithm 1 · p. 9
  • FFT and inverse FFT are typically the most costly operations; their cost scales as n log(n), but efficiency depends strongly on grid-point count.Slightly changing the grid size can sometimes yield orders-of-magnitude speed-ups; powers of two are optimal, though not always feasible.Source: Algorithm 1 · p. 9
  • The constitutive model is evaluated at nodes using the nodal deformation gradient, so FEM integration-point constitutive implementations can also be used.Source: Algorithm 1 · p. 9

4 Application to Hyper-elasticity, formulated in the reference configuration

The paper demonstrates its finite-strain implementation on a hyper-elastic two-phase cell under large simple shear. The example resolves phase-wise stress and deformation, while the authors delimit the model and grid settings in which the results apply.

4.2 Constitutive model

  • The reference-configuration hyper-elastic model uses a linear relation between second Piola-Kirchhoff stress S and Green-Lagrange strain E through the elastic stiffness tensor C.C is the standard fourth-order elastic stiffness tensor, parameterized using Lamé constants.Source: §4.2 · p. 9
  • The model is not a proper elastic constitutive model for very large deformations, and is used as a simple, sufficiently realistic approximation within the analyzed deformation range.The authors report perfect agreement with Kabel et al.'s analytical 2-D laminate solution, though that comparison is not shown.Source: §4.2 · p. 9
  • The reference-configuration constitutive law is connected to the deformation-gradient formulation through definitions of Green-Lagrange strain and first Piola-Kirchhoff stress.Source: §4.2 · p. 10

4.3 Consistent tangent

  • The consistent tangent is derived by linearizing the stress relation and Green-Lagrange strain, then combining those linearizations and using the stiffness tensor's right symmetry.Source: §4.3 · p. 10

4.4 Example

  • The example considers a three-dimensional two-phase unit cell with a cubic particle whose stiffness is 10 times higher than the matrix, loaded in simple shear with macroscopic γ̄ = 1.The parameters include ν = 0.3.Source: §4.4 · p. 10
  • The hard phase carries significantly higher equivalent second Piola-Kirchhoff stress, while the soft phase accommodates most deformation, with no nodal artifacts despite discontinuous deformation gradients across the interface.The unit-cell example runs in under a minute on an ordinary desktop and requires five iterations in one increment.Source: §4.4 · p. 10
  • The regular grid required by FFT can constrain problem choice; for inclusion volume-fraction homogenization, coarse and moderately fine results may be grid-sensitive but converge with refinement.Source: §4.4 · p. 10

5 Application to Simo elasto-plasticity, defined in the current configuration

The section applies finite-strain FFT methods to Simo elasto-plasticity in the current configuration, then demonstrates the implementation on a dual-phase steel micrograph. It reports localized plastic strain and contrast-dependent solver runtime.

5.1 Introduction

  • The Simo model demonstrates that a current-configuration material model can be treated using finite-element procedures, with its implementation remaining below 200 lines of Python.The code is longer than the hyperelastic example because its stress update and constitutive tangent are more involved.Source: §5.1 · p. 10

5.2 Constitutive model

  • The elasto-plastic Simo model uses a multiplicative elastic-plastic split, logarithmic elastic strain, and implicit plasticity with an elastic-predictor plastic-corrector algorithm [30].Its yield stress hardens with accumulated plastic strain, and plastic flow follows normality.Source: §5.2 · p. 11

5.3 Consistent tangent

  • The consistent tangent combines geometric nonlinearity with the linearized constitutive response, including the strain definition; its discrete elasto-plastic switch affects convergence around yield.An explicit expression for the constitutive tangent is available in references.Source: §5.3 · p. 11

5.4 Pull-back to the undeformed configuration

  • The pull-back maps the Kirchhoff stress and consistent tangent to the undeformed configuration without introducing additional approximations.The tangent is defined in the last-known iterative configuration.Source: §5.4 · p. 12

5.5 Example

  • The steel example uses a thresholded 451×451-pixel micrograph under plane strain, with martensite and ferrite treated as elastically homogeneous phases.The example uses linear hardening and χ = 2.Source: §5.5 · p. 12
  • At macroscopic equivalent strain 0.2, plastic strain peaks at 0.51 near hard phases; the 250-increment run averages 2.2 Newton iterations per increment.The plastic strain is higher in the soft phase, with localization aligned at ±45° to the applied shear.Source: §5.5 · p. 12, Figure 3 · p. 13
  • For phase contrasts χ = 2, 4, and 8, Newton iteration counts are practically independent of contrast, while runtime is lower at lower χ.The authors attribute the runtime trend to more homogeneous constitutive K entries at lower contrast; the CG iteration count is unavailable from the standard Python routine.Source: §5.5 · p. 12, Figure 4 · p. 13

6 Summary

The method solves large-deformation equilibrium in periodic cells using a variational formulation, globally supported trigonometric interpolation, and a projection operator for compatibility. The authors report a general, memory-efficient approach whose FFT-based computational efficiency and compact Python implementation support use with different constitutive models.

6 Summary

  • The method solves large-deformation equilibrium in periodic cells through a variational formulation, trigonometric interpolation, and projection-based compatibility enforcement.It uses nodal deformation gradients and first Piola-Kirchhoff stresses rather than nodal displacements and forces.Source: §6 · p. 14
  • The method is suitable for arbitrary constitutive models in a finite-deformation setting.The constitutive implementation is independent of the solution strategy and can be used with FFT or conventional finite-element solvers.Source: §6 · p. 14
  • The method is memory efficient and avoids involved data storage.Source: §6 · p. 14
  • Computational efficiency results from using the Fast Fourier Transform.Source: §6 · p. 14
  • A 59-line Python code demonstrates the method for a three-dimensional finite-strain problem.Source: §6 · p. 14

A Nomenclature

The nomenclature specifies tensor orders, contraction notation, identity tensors, and a sequence notation used in the paper.

A Nomenclature

  • The nomenclature identifies A_ij as a second-order tensor and A_ijkl as a fourth-order tensor.Source: A Nomenclature · p. 15
  • It gives double-dot contraction notation for second- and fourth-order tensors.The listed forms are C = A : B, with component expressions for fourth-order–second-order and fourth-order–fourth-order contractions.Source: A Nomenclature · p. 15
  • It names I_RT as the fourth-order right-transposed identity, 2(I + I_RT) as the fourth-order symmetrization tensor, and A_1, A_2, …, A_n as a sequence.Source: A Nomenclature · p. 15

B Python code

The Python appendix presents code for a three-dimensional example, using tensor matrices and standard scientific Python modules. It also introduces the problem dimensions as the spatial dimension count and grid-point count in each dimension.

B Python code

  • The appendix code stores quantities as three-dimensional matrices of second- or fourth-order tensors, avoiding sparse data structures and storing only what is needed.Defining the Fourier coefficients of the projection operator is identified as the somewhat costly part; a vectorized version is available online.Source: B Python code · p. 16

Import modules

  • The import section uses NumPy and SciPy’s sparse linear algebra module, with NumPy providing linear algebra operations for scientific computation.Source: Import modules · p. 16

Problem dimensions

  • The problem dimensions are the number of spatial dimensions and the number of grid points in each spatial dimension.Source: Problem dimensions · p. 16

Tensor products

The code defines tensor operations with NumPy’s index notation and initializes identity tensors used to evaluate the constitutive response.

Tensor products

  • The tensor-product section uses NumPy einsum index notation to define tensor operations and products.The listed operations include transposition, contractions, products, and dyadic products.Source: Tensor products · p. 16

Define tensors

  • The code initializes second- and fourth-order identity tensors, including a symmetrized fourth-order identity, for constitutive-response evaluation.Source: Define tensors · p. 17

Projection

The implementation constructs Fourier-space projection operations, defines a spatially varying constitutive model, and solves for equilibrium with Newton iterations and conjugate gradients.

Projection

  • Projection operations apply Fourier transforms, contract fields with Fourier coefficients, and inverse-transform the result; the coefficients are initialized separately.The implementation also uses fftshift and its inverse to align the zero frequency with the grid center.Source: Projection · p. 17

Problem definition and constitutive model

  • The example defines bulk- and shear-modulus grids for a cubic particle in a homogeneous matrix and evaluates stress and tangent tensors from the deformation gradient.The constitutive function evaluates first Piola-Kirchhoff stress P and tangent stiffness K across the grid.Source: Problem definition and constitutive… · p. 18

Newton iterations

  • Starting from P = O and F = I, the solver applies the macroscopic deformation, then uses Newton iterations and conjugate gradients to update the field toward equilibrium.Each iteration evaluates P and K to form the residual and tangent before solving for the deformation-gradient update.Source: Newton iterations · p. 18

C Even-sized grid

For even-sized grids, Nyquist frequencies prevent compatibility and stress equilibrium from being satisfied simultaneously; the implementation sets the projection to zero at those frequencies.

C Even-sized grid

  • Even-sized grids cannot satisfy both deformation-gradient compatibility and stress equilibrium simultaneously because of Nyquist frequencies.Setting the projection to zero at Nyquist frequencies recovers compatibility; setting it to identity recovers equilibrated stress.Source: C Even-sized grid · p. 20
  • The implementation sets the projection to zero at Nyquist frequencies, accepting approximate equilibrium when exact equilibrium is not generally a concern.With this choice, comparison against Kabel et al. produced perfect agreement with the analytical solution.Source: C Even-sized grid · p. 20
  • For the Section 5 example, results using an even grid with one additional pixel on each side were identical to those from the other approach.Source: C Even-sized grid · p. 20

D Projection operator

The appendix explains the projection operator in Fourier space and establishes that it maps tensor fields onto the zero-mean compatible fields.

D Projection operator

  • In Fourier space, the projection acts row-wise, projecting each row of the input onto the direction of the scaled frequency.The argument checks that the projected rows have vanishing curl, making them gradients of scalar potentials.Source: D Projection operator · pp. 20–21
  • The appendix also establishes that the operator's range includes the whole set of compatible fields, representing them as gradients of scalar potentials.Source: D Projection operator · p. 21
  • Operator G maps square-integrable second-order tensor fields onto the zero-mean compatible fields.Source: D Projection operator · p. 21
Loading 1603.08893v3…