FFT-based Solver Research Articles

Mérope: A microstructure generator for simulation of heterogeneous materials

Mérope is a software devoted to the geometrical design and the discretization of microstructures of random heterogeneous materials. Mérope aims at building large samples of microstructured materials, called Representative Volume Elements, in order to derive their effective physical behaviors. Various microstructures are supported: spherical, polyhedral or spheropolyhedral inclusions, polycrystals, Gaussian fields and Boolean combinations of these. Discretization takes two forms: either regular Cartesian grids of (composite) voxels for computations with FFT-based solvers, or tetrahedral meshes for computations with Finite Element solvers. A special emphasis on the code has been put on performance, which will be further improved in the future.This article aims at introducing the main features of the software as well as exemplifying its use.

Accelerating the force-coupling method for hydrodynamic interactions in periodic domains

The efficient simulation of fluid-structure interactions at zero Reynolds number requires the use of fast summation techniques in order to rapidly compute the long-ranged hydrodynamic interactions between the structures. One approach for periodic domains involves utilising a compact or exponentially decaying kernel function to spread the force on the structure to a regular grid where the resulting flow and interactions can be computed efficiently using an FFT-based solver. A limitation to this approach is that the grid spacing must be chosen to resolve the kernel and thus, these methods can become inefficient when the separation between the structures is large compared to the kernel width. In this paper, we address this issue for the force-coupling method (FCM) by introducing a modified kernel that can be resolved on a much coarser grid, and subsequently correcting the resulting interactions in a pairwise fashion. The modified kernel is constructed to ensure rapid convergence to the exact hydrodynamic interactions and a positive-splitting of the associated mobility matrix. We provide a detailed computational study of the methodology and establish the optimal choice of the modified kernel width, which we show plays a similar role to the splitting parameter in Ewald summation. Finally, we perform example simulations of rod sedimentation and active filament coordination to demonstrate the performance of fast FCM in application.

Implementation of a phase field damage model with a nonlinear evolution equation in an FFT-based solver

This paper focuses on the numerical implementation of phase-field models of fracture using the Fast Fourier Transform (FFT) based numerical method. Recent studies on phase-field models focus on the discussions of the choice of the value of regularization length proposed to smear the discontinuity of the sharp crack. Some studies argue that it should be considered as a material property because it has a significant impact on the mechanical behavior of a material in some phase-field models, for instance, in the model proposed by Miehe. However, our results in this study for heterogeneous materials have shown that the choice of regularization length not only affects the macroscopic mechanical behavior but also the local crack propagation patterns. As a result, it can be challenging to select an appropriate value that produces both accurate macroscopic responses and local crack patterns for certain phase-field models, such as Miehe’s model. Thus, the phase-field model proposed by Wu, which has been proven to reduce the length sensitivity for homogeneous material, has been successfully implemented in an FFT-based solver with the application of Newton–Krylov algorithm.The length sensitivity of Wu’s phase-field model for heterogeneous materials is also investigated in this paper. Our tests show that Wu’s phase-field model has partial length sensitivity for heterogeneous materials. However, the main reason for this sensitivity is that one phase enters the damage zone of other phase, which is different from Miehe’s model. As a result, a set of criteria for accurately choosing the regularization value has been established for Wu’s model in this work. Meanwhile, it has also been found that Wu’s model may be more suitable than Miehe’s model for brittle failure due to the introduction of an elastic stage.

Accounting for weak interfaces in computing the effective crack energy of heterogeneous materials using the composite voxel technique

AbstractWe establish a computational methodology to incorporate interfaces with lower crack energy than the surrounding phases when computing the effective crack energy of brittle composite materials. Recent homogenization results for free discontinuity problems are directly applicable to the time-discretized Francfort-Marigo model of brittle fracture in the anti-plane shear case, and computational tools were introduced to evaluate the effective crack energy on complex microstructures using FFT-based solvers and a discretization scheme based on a combinatorially consistent grid. However, this approach only accounts for the crack resistance per volume and is insensitive to the crack resistance of the interface which is expected to play a significant role by considerations from materials science. In this work we introduce a remedy exploiting laminate composite voxels. The latter were originally introduced to enhance the accuracy of solutions for elasticity problems on regular voxel grids. We propose an accurate approximation of the effective crack energy of a laminate with weak interface where an explicit solution is available. We incorporate this insight into an efficient algorithmic framework. Finally, we demonstrate the capabilities of our approach on complex microstructures with weak interfaces between different constituents.

A Fast Hybrid Pressure-Correction Algorithm for Simulating Incompressible Flows by Projection Methods

To enforce the conservation of mass principle, a pressure Poisson equation arises in the numerical solution of incompressible fluid flow using the pressure-based segregated algorithms such as projection methods. For unsteady flows, the pressure Poisson equation is solved at each time step usually in physical space using iterative solvers, and the resulting pressure gradient is then applied to make the velocity field divergence-free. It is generally accepted that this pressure-correction stage is the most time-consuming part of the flow solver and any meaningful acceleration would contribute significantly to the overall computational efficiency. The objective of the present work was to develop a fast hybrid pressure-correction algorithm for numerical simulation of incompressible flows around obstacles in the context of projection methods. The key idea is to adopt different numerical methods/discretisations in the sub-steps of projection methods. Here, a classical second-order time-marching projection method, which consists of two sub-steps, was chosen for the purposes of demonstration. In the first sub-step, the momentum equations were discretised on unstructured grids and solved by conventional numerical methods, here a meshless method. In the second sub-step (pressure-correction), the proposed algorithm adopts a double-discretisation system and combines the weighted least-squares approximation with the essence of immersed boundary methods. Such a design allowed us to develop an FFT-based solver to speed up the solution of the pressure Poisson equation for flow cases with obstacles, while keeping the implementation of the boundary conditions for the momentum equations as easy as conventional numerical methods do with unstructured grids. The numerical experiments of five test cases were performed to verify and validate the proposed hybrid algorithm and evaluate its computational performance. The results showed that the new FFT-based hybrid algorithm works and is robust, and it was significantly faster than the multigrid-based reference method. The hybrid algorithm opens an avenue for the development of next-generation high-performance parallel computational fluid dynamics solvers for incompressible flows.

Investigations on the influence of the boundary conditions when computing the effective crack energy of random heterogeneous materials using fast marching methods

Recent stochastic homogenization results for the Francfort–Marigo model of brittle fracture under anti-plane shear indicate the existence of a representative volume element. This homogenization result includes a cell formula which relies on Dirichlet boundary conditions. For other material classes, the boundary conditions do not effect the effective properties upon the infinite volume limit but may have a strong influence on the necessary size of the computational domain. We investigate the influence of the boundary conditions on the effective crack energy evaluated on microstructure cells of finite size. For periodic boundary conditions recent computational methods based on FFT-based solvers exploiting the minimum cut/maximum flow duality are available. In this work, we provide a different approach based on fast marching algorithms which enables a liberal choice of the boundary conditions in the 2D case. We conduct representative volume element studies for two-dimensional fiber reinforced composite structures with tough inclusions, comparing Dirichlet with periodic boundary conditions.

FFT-based investigation of the shear stress distribution in face-centered cubic polycrystals

The onset of nonlinear effects in metals, such as plasticity and damage, is strongly influenced by the heterogeneous stress distribution at the grain level. This work is devoted to studying the local stress distribution in fcc polycrystals using FFT-based solvers. In particular, we focus on the distribution of shear stresses resolved in the slip systems as the critical driving force for plastic deformations. Specific grain orientations with respect to load direction are investigated in the linear elastic regime and at incipient plastic deformations based on a large ensemble of microstructures. The elastic anisotropy of the single crystal is found to have a crucial influence on the scatter of the stress distribution, whereas the Young’s modulus in the respective crystal direction governs the mean stress in the grain. It is further demonstrated that, for higher anisotropy, the shear stresses deviate from the normal distribution and are better approximated by a log-normal fit. Comparing the full-field simulations to the Maximum Entropy Method (MEM), reveals that the MEM provides an excellent prediction up to the second statistical moment in the linear elastic range. In a study on the spatial distribution of shear stresses, the grain boundary is identified as a region of pronounced stress fluctuations and as the starting point of yielding during the elastic–plastic transition.

A simple extension of FFT-based methods to strain gradient loadings - Application to the homogenization of beams and plates with linear and non-linear behaviors

Because of their simplicity, efficiency and ability for parallelism, FFT-based methods are very attractive in the context of numerical periodic homogenization, especially when compared to standard FE codes used in the same context. They allow applying to a unit-cell a uniform average strain with a periodic strain fluctuation that is an unknown quantity. Solving the problem allows to evaluate the complete stress-strain fields. The present work proposes to extend the use of the method from uniform loadings (i.e. uniform applied strain) to strain gradient loadings (i.e. strain fields with a uniform strain gradient) while keeping the algorithm as simple as possible. The identification of a subset of strain gradient loadings allows for a minimally invasive modification of the iterative algorithm so that the implementation is straightforward. In spite of a reduced subset of 9 independent loadings among the 18 available, the second part of the paper demonstrates that it is enough for considering the homogenization of beams and plates. A first application validates the approach and compares it to another FFT-based method dedicated to the homogenization of plates. The second application concerns the homogenization of beams, for the first time considered (to author's knowledge) with an FFT-based solver. The method applies to different beam cross-sections and the proposition of using composite voxels drastically improves the numerical solution when the beam cross-section is not conform with the spatial discretization, especially for torsion loading. As a result, the massively parallel AMITEX_FFTP code has been slightly modified and now offers a new functionality, allowing the users to prescribe torsions and flexions to beam or plate heterogeneous unit-cells.

A dual resolution phase-field solver for wetting of viscoelastic droplets.

We present a new and efficient phase‐field solver for viscoelastic fluids with moving contact line based on a dual‐resolution strategy. The interface between two immiscible fluids is tracked by using the Cahn‐Hilliard phase‐field model, and the viscoelasticity incorporated into the phase‐field framework. The main challenge of this approach is to have enough resolution at the interface to approach the sharp‐interface methods. The method presented here addresses this problem by solving the phase field variable on a mesh twice as fine as that used for the velocities, pressure, and polymer‐stress constitutive equations. The method is based on second‐order finite differences for the discretization of the fully coupled Navier–Stokes, polymeric constitutive, and Cahn–Hilliard equations, and it is implemented in a 2D pencil‐like domain decomposition to benefit from existing highly scalable parallel algorithms. An FFT‐based solver is used for the Helmholtz and Poisson equations with different global sizes. A splitting method is used to impose the dynamic contact angle boundary conditions in the case of large density and viscosity ratios. The implementation is validated against experimental data and previous numerical studies in 2D and 3D. The results indicate that the dual‐resolution approach produces nearly identical results while saving computational time for both Newtonian and viscoelastic flows in 3D.

On the impact of the mesostructure on the creep response of cellular NiAl-Mo eutectics

Directionally solidified eutectics of NiAl matrix and fibrous refractory metals, like Mo, can form cellular mesostructures with significant fiber misalignment and changing fiber volume fraction, for example, when being solidified at high growth rates or when increased solidification intervals are present in the alloys. In order to reveal the deteriorating impact of the mesostructure, i.e., the volume fraction and aspect ratio of the well-aligned cells, on the creep response of such cellular eutectics, we rely on scale-bridging numerical simulations, using the level-set framework by Sonon et al. [1] for microstructure generation and FFT-based solvers for computing the creep response. Our results indicate, firstly, that the fraction of properly aligned regions in cellular NiAl composites is lower than estimated in earlier experimental studies, due to the existence of degenerated regions surrounding the well-aligned cell interiors. Secondly, studying the influence of the cell aspect ratio shows that the apparent stress exponent of the composite is very sensitive with respect to this parameter, providing a possible explanation for the large scatter of experimentally determined stress exponents in previous studies. A comparison of the numerical simulations to a linear rule of mixtures and the frequently applied analytical Kelly-Street model illustrates that both fail to accurately describe the magnitude of minimum creep rates in the investigated ranges of volume fractions and aspect ratios. The heterogeneity of the strain-rate field on the mesoscale is identified as the primary error source, demonstrating that either numerical simulations or more sophisticated analytical models are required for reliably predicting for the creep response of cellular materials.

Analysis of a model of field crack mechanics for brittle materials

A computational model for arbitrary brittle crack propagation, in a fault-like layer within a 3-d elastic domain, and its associated quasi-static and dynamic fields is developed and analyzed. It uses a FFT-based solver for the balance of linear momentum and a Godunov-type projection-evolution method for the crack evolution equation. As applications, we explore the questions of equilibria and irreversibility for crack propagation with and without surface energy, existence of strength and toughness criteria, crack propagation under quasi-static and dynamic conditions, including Modes I, II and III, as well as multiaxial compressive loadings.

An FFT solver used for virtual Dynamic Mechanical Analysis experiments: Application to a glassy/amorphous system and to a particulate composite

FFT-based solvers are increasingly used by many researcher groups interested in modelling the mechanical behavior associated to a heterogeneous microstructure. A development is reported here that concerns the viscoelastic behavior of composite structures generally studied experimentally through Dynamic Mechanical Analysis (DMA). A parallelized computation code developed with complex-valued quantities provides virtual DMA experiments directly in the frequency domain on a heterogeneous system described by a voxel grid of mechanical properties. The achieved precision and computation times are very good. An effort has been made to illustrate the application of such a virtual DMA tool through two examples from the literature: the modelling of glassy/amorphous systems at a small scale and the modelling of experimental data obtained in temperature sweeping mode by DMA on a particulate composite made of glass beads and a polystyrene matrix, at a larger scale. Both examples show how virtual DMA can contribute to question, analyze, and understand relaxation phenomena on either theoretical or experimental points of view.

FFT-based simulations of slip and kink bands formation in 3D polycrystals: Influence of strain gradient crystal plasticity

The ability of softening strain gradient plasticity to simulate intragranular plastic slip localization modes called slip and kink bands, at is investigated at incipient plasticity. A strain gradient crystal plasticity model based on a quadratic energy contribution of the Nye tensor has been implemented within the massively parallel FFT-based solver AMITEX_FFTP. It is shown that a classical evaluation of the double curl operator does not properly account for the backstress induced by GNDs in slip and kink bands. Consequently, a better suited differentiation operator is proposed for this purpose. Besides, four practical implementations of interface conditions at grain boundaries are implemented and discussed. It is demonstrated that these conditions have a strong influence on GND induced hardening but weakly affect the formation of slip localization patterns.A theoretical and numerical study of an idealized kink band shows that this modeling framework, contrary to classical crystal plasticity models, is able to capture their essential physical features. A quantitative analysis of high resolution 3D polycrystalline simulations, based on imaged processing and a peak detection on plastic slip profiles, is applied to show that gradient effects strongly influence the mechanisms of selection between slip and kink banding, to the favor of the first. In particular, some kink bands are shown to decompose into an homogeneous distribution of parallel slip bands of finite length. Finally, the competition between structural effects and energetic cost of lattice curvature in the formation of localized slip bands networks in polycrystals is discussed. The present results demonstrate the relevance of softening crystal plasticity modeling including an energetic contribution of GND density to accurately capture slip localization in polycrystals.

A modified FFT-based solver for the mechanical simulation of heterogeneous materials with Dirichlet boundary conditions

A modified FFT-based solver for the mechanical simulation of heterogeneous materials with Dirichlet boundary conditions

Periodic smoothing splines for FFT-based solvers

The aim of this paper is to develop a periodic smoother based on splines for FFT-based solvers. Spurious oscillations in FFT-based methods are shown to be due to pseudo-spectral differentiation of discontinuous fields. An automatic smoother based on polynomial splines is developed, which permits to add smoothness to initial material properties. The method, which is applied in various problems including conductivity, elasticity and field dislocation mechanics, improves significantly the local fields and reduces spurious oscillations.

A micromechanical analysis of intergranular stress corrosion cracking of an irradiated austenitic stainless steel

Irradiation Assisted Stress Corrosion Cracking (IASCC) is a material degradation phenomenon affecting austenitic stainless steels used in nuclear Pressurized Water Reactors (PWR), leading to the initiation and propagation of intergranular cracks. Such phenomenon belongs to the broader class of InterGranular Stress Corrosion Cracking (IGSCC). A micromechanical analysis of IGSCC of an irradiated austenitic stainless steel is performed in this study to assess local cracking conditions. A 304L proton irradiated sample tested in PWR environment and showing intergranular cracking is investigated. Serial sectioning, Electron BackScatter Diffraction (EBSD) and a two-step misalignment procedure are performed to reconstruct the 3D microstructure over an extended volume, to assess statistically cracking criteria. A methodology is also developed to compute Grain Boundary (GB) normal orientations based on the EBSD measurements. The statistical analysis shows that cracking occurs preferentially for GB normals aligned with the mechanical loading axis, but also for low values of the Luster-Morris slip transmission parameter. Micromechanical simulations based on the reconstructed 3D microstructure, FFT-based solver and crystal plasticity constitutive equations modified to account for slip transmission at grain boundaries are finally performed. These simulations rationalize the correlation obtained experimentally into a single stress-based criterion. The actual strengths and weaknesses of such micromechanical approach are discussed.

Improving the initial guess for the Newton-Raphson protocol in time-dependent simulations

A general linearisation procedure for the consistent tangent of a small-strain visco-plastic material model is presented in this note. The procedure is based on multi-variable linearisation around a so-called 'reference state'. In particular, the linerarisation of the time integration scheme is found to yield an extra term compared to classical expressions, which only appears because the material response is time-dependent. It has the effect of yielding a very accurate initial guess for the Newton-Raphson protocol based on the ongoing viscous flow. It is shown, using a modern variational FFT-based solver, that the extra term reduces both the CPU time and the number of Newton-Raphson iterations by around a factor two.

Non-linear composite voxels for FFT-based explicit modeling of slip bands: Application to basal channeling in irradiated Zr alloys

Explicitly modeling slip bands is a promising method to simulate the consequences of intragranular plastic slip localization on the behavior of strongly localizing crystals such as irradiated metals. In this study, we propose a very efficient framework to implement this modeling strategy based on a massively parallel FFT-based solver enhanced with composite voxels to simulate polycrystals with a regular lamellar structure alternating plastic slip bands embedded within larger elastic bands. To this end, generic composite voxel models have been formulated and implemented, allowing to extend the approach to any type of constitutive behavior in the future. In particular, we show that laminate composite voxels are particularly well suited to improve standard FFT-based solvers in this context, and enable to run polycrystalline simulations accurately accounting for slip bands whose width is two orders of magnitude below the grain size. To demonstrate the potential of this approach, we apply it to an idealized material designed to approximate basal channeling in irradiated textured Zr alloys. Varying the width and the spacing of the modeled slip bands, we show that it allows to capture the increase of grain boundary normal stress concentrations or the enhanced Bauschinger effect induced by an increased localization of plastic slip.

An immersed boundary projection method for simulating the inextensible vesicle dynamics

We develop an immersed boundary projection method (IBPM) based on an unconditionally energy stable scheme to simulate the vesicle dynamics in a viscous fluid. Utilizing the block LU decomposition of the algebraic system, a novel fractional step algorithm is introduced by decoupling all solution variables, including the fluid velocity, fluid pressure, and the elastic tension. In contrast to previous works, the present method preserves both the fluid incompressibility and the interface inextensibility at a discrete level simultaneously. In conjunction with an implicit discretization of the bending force, the present method alleviates the time-step restriction, so the numerical stability is assured by non-increasing total discrete energy during the simulation. The numerical algorithm takes a linearithmic complexity by using preconditioned GMRES and FFT-based solvers. The grid convergence studies confirm the solution variables exhibit first-order convergence rate in L2-norm. We demonstrate the numerical results of the vesicle dynamics in a quiescent fluid, Poiseuille flow, and shear flow, which are congruent with the results in the literature.

FFT-based solver for higher-order and multi-phase-field fracture models applied to strongly anisotropic brittle materials

This paper presents the application of a fast Fourier transform (FFT) based method to solve two phase field models designed to simulate crack growth of strongly anisotropic materials in the brittle regime. By leveraging the ability of the FFT-based solver to generate solutions with higher-order and global continuities, we design two simple algorithms to capture the complex fracture patterns (e.g. sawtooth, and curved crack growth) common in materials with strongly anisotropic surface energy via the multi-phase-field and high-order phase-field frameworks. A staggered operator-split solver is used where both the balance of linear momentum and the phase field governing equations are formulated in the periodic domain. The unit phase field of the initial failure region is prescribed by the penalty method to alleviate the sharp material contrast between the initial failure region and the base material. The discrete frequency vectors are generalized to estimate the second and fourth order gradients such that the Gibbs effect near shape interfaces or jump conditions can be suppressed. Furthermore, a preconditioner is adopted to improve the convergence rate of the iterative linear solver. Three numerical experiments are used to systematically compare the performance of the FFT-based method in the multi-phase-field and high-order phase-field frameworks.