Phase field modelling of fracture and fatigue in Shape Memory Alloys

To constitutively describe the material behaviour of SMAs we follow the so-called unified model by Lagoudas and co-workers [39, 40, 1], including recent extensions to capture gradual phase transformations, and the stress-dependencies of the inelastic recoverable strain and the phase diagram slope [41, 42]. The model builds upon the rule of mixtures to determine the magnitude of relevant properties for material points located in the phase transformation region. Thus, the effective value of any relevant phase-dependent parameter ($\Theta$) is a function of the martensitic volume fraction $\xi$ ($0\leq\xi\leq 1$) and its magnitude in the austenitic ($\Theta^{A}$) and martensitic ($\Theta^{M}$) phases;

Fracture and fatigue in SMAs is predicted by using a variational phase field model [19, 43, 44]. Since the early work by Francfort and Marigo [18], phase field fracture models aim at providing a variational framework for the concept of crack advance driven by the competition between toughness (surface energy density, $G_{c}$) and energy release rate $G$; as first proposed by Griffith [17] for elastic solids, and later extended to account for inelastic energy dissipation by Orowan [45]. In the present study, we will consider both a constant material toughness $G_{c}$ and the case in which the critical energy release rate is determined from its austenite and martensite counterparts, $G_{c}(\xi)$, using the rule of mixtures.

The discrete crack is approximated through an auxiliary field variable $\phi$, which varies between $\phi=0$, intact material, and $\phi=1$, fully cracked material. The size of the regularized crack surface is governed by the choice of $\ell$, the phase field model-inherent length scale. Thus, the fracture energy due to the creation of a crack can be approximated as:

The main difference between them is the fact that the AT1 model has a non-zero elastic limit. Other constitutive choices have been proposed, such as the so-called phase-field regularized cohesive zone model (PF-CZM) [46]; the reader is referred to Ref. [47] for a detailed numerical comparison. To retain generality, we proceed to present our phase field formulation without making specific constitutive choices for $w\left(\phi\right)$ and $c_{w}$, as both the AT1 and the AT2 models will be considered in this study.

The total potential energy can be expressed as a function of the contributions from the mechanical, thermal and fracture terms as:

Consider now the total potential energy of the solid, Eq. (19). The strong form can be readily derived by taking the first variation with respect to $\bm{\varepsilon}$, $T$ and $\phi$, and making use of Gauss’ divergence theorem. Thus, the coupled field equations read,

The implementation of the constitutive model follows the work by Qidwai and Lagoudas [52], where the evolution of the transformation strain tensor is incrementally integrated using the backward Euler method. Thus, we use the common notation of adding the subscript $n+1$ to quantities in the current time step while the subscript $n$ denotes the previous time step, e.g. for the transformation strain tensor:

In addition, as the problem is solved in an iterative manner, the superscript $(k)$ is used for the iteration counter; i.e., (26) becomes

The incremental stresses are computed using a return mapping algorithm; a purely elastic trial state is followed by a transformation correction stage. Thus, assuming a uniform temperature and denoting $\mathcal{\bm{C}}$ as the elastic stiffness matrix, the elastic stress prediction is given by

The return mapping algorithm is then used to enforce satisfying the transformation consistency condition $\dot{\Phi}=0$. Specifically, the convex cutting plane algorithm [53] is used, which differs from the commonly used closest point projection algorithm as follows. Considering the flow rule (5) and (27), the incremental transformation tensor reads,

In the convex cutting plane algorithm the implicit dependence on the transformation direction $\bm{\Lambda}^{t}\left(\bm{\sigma}_{n+1}^{(k+1)}\right)$ is relaxed, and (29) can be reformulated as,

The total current strain is held constant during the iterative correction such that, considering (12), the incremental Cauchy stresses read,

In an iterative setting, the consistency condition implies,

And finally, solving (33) for the correction in the martensitic volume fraction at iteration $(k)$ gives,

The transformation strain tensor can then be updated via (27) and (30). The updated values of the transformation strain and elastic stiffness are used to calculate an updated stress tensor via,

The last step involves computing the consistent material Jacobian $\mathcal{\bm{L}}$ [54, 52, 42], which for isothermal conditions is given by,

Fracture is assumed to be driven by the total strain energy density $\psi$; i.e., both the elastic $\psi^{e}$ and transformation $\psi^{t}$ strain energy densities contribute to cracking on equal footing. This choice is mainly phenomenological but has some physical background, as it appears sensible to assume a contribution from the transformation strains to the fracture process. Similar approaches have been adopted by other authors in relation to other inelastic quantities; see e.g., Ref. [21, 55] for examples in the context of plasticity and ductile damage. To maintain resistance in compression and during crack closure, the elastic contribution to the strain energy density is decomposed into volumetric and deviatoric parts, following Amor et al. [56]. In this approach, the deviatoric and tensile volumetric components contribute to fracture but the compressive volumetric term does not. Thus, the elastic strain energy is decomposed into the following two terms:

We proceed to formulate the two-field weak form of the problem and subsequently derive the stiffness matrices and residuals applying a finite element discretisation. Thus, consider the total potential energy of the solid, Eq. (19), under isothermal conditions and in the absence of body forces and external tractions. The first variation of (19) with respect to $\bm{\varepsilon}$ and $\phi$, gives

Now make use of Voigt notation and consider a 3D solid. The displacement field $\bm{u}$ and the phase field $\phi$ are discretised as

We proceed to formulate the residuals and the stiffness matrices considering this finite element discretisation and the fact that (42)-(43) must hold for arbitrary values of $\delta\bm{u}$ and $\delta\phi$. The associated discrete equations can be formulated as the following residuals with respect to the displacement field and the phase field variable, respectively,

A global iterative scheme is adopted to obtain the displacement $\bm{u}$ and phase field $\phi$ solutions for which $\bm{r}^{\bm{u}}=\bm{0}$ and $\textbf{r}^{\phi}=\bm{0}$;

We develop a numerical implementation that can accommodate both staggered [43] and monolithic quasi-Newton schemes [58, 23]. As shown by Wu et al. [58] for the PF-CZM model and by Kristensen and Martínez-Pañeda [23] for the AT2 model, the use of quasi-Newton methods such as the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm enables a robust implementation of monolithic schemes that retain unconditional stability, speeding up calculations by several orders of magnitude. In this work, this is particularly relevant for the analysis of fatigue, as otherwise calculations would be prohibitive [23]. The reader is referred to Refs. [43] and [23] for details on the implementation of staggered and monolithic quasi-Newton solution schemes, respectively.

Thus, the nodal displacements in the outer boundary of the finite element model can be prescribed to evaluate the crack tip behaviour at a given value of $K_{I}$. However, one should note that this is under the assumption that the inelastic region is small, generally referred to as small transformation zone conditions; analogous to the so-called small scale yielding conditions in elastic-plastic materials. Accordingly, Eqs. (51)-(52) make use of the elastic constants for the austenitic phase.

We proceed to model the failure of a cracked square plate subjected to tension, a paradigmatic benchmark in phase field fracture [43, 57, 29, 23]. The specimen has an initial horizontal crack going from the left side to the center of the specimen, both vertical and horizontal displacements are restricted in the bottom boundary, and we load the plate by prescribing the vertical displacement in the upper edge, $u^{\infty}$. The geometric setup, dimensions (in mm) and boundary conditions are given in Fig. 6. The constitutive behaviour of the material is characterised by the parameters listed in Table 1 while the phase field model uses $G_{c}=4.1$ kJ/m² and $\ell=0.0075$ mm. The characteristic element size is at least 7 times smaller than $\ell$ along the extended crack plane. A total of 45,571 quadrilateral linear elements with full integration are used.

The results obtained are shown in Fig. 7 in terms of the force versus displacement response, with contours of martensite volume fraction $\xi$ and phase field fracture parameter $\phi$ at different stages embedded in the figure. Differences from the response commonly observed in this classic benchmark are notable. In a linear elastic homogeneous material the plate fails in an unstable manner, with the force versus displacement curve exhibiting a very large drop immediately after reaching the peak load - see, e.g. [43, 23]. Contrarily, in the SMA sample a significant toughening effect is observed; there is an initial drop in the load associated with the first instance of crack growth but then the crack progresses in a stable manner until complete failure of the sample. The toughening effect observed is undoubtedly related to the energy dissipated due to transformation. As shown in Fig. 7, a large transformation zone develops in the sample, exhibiting a shear banding-like behaviour that resembles that observed in elastic-plastic materials [68]. However, one should note that, for the material properties here considered, the crack still follows the mode I fracture path.

An asymmetrically notched plane strain bar is investigated to model mixed-mode fracture and crack coalescence. The double-notched bar, depicted in Fig. 8, is clamped at the bottom end ($u_{x}=u_{y}=0$) and subjected to a vertical displacement $u^{\infty}$ at the top edge. Two circular notches of radii 2.5 mm have been geometrically introduced. The bar is assumed to be made of the equiatomic nitinol SMA whose properties and model parameters are listed in Table 1. The phase field length scale equals 0.2 mm and the finite element mesh is chosen accordingly, with the characteristic element length in the region between the two notches being on the order of 0.05 mm. A total of 23,622 quadrilateral plane strain linear elements are used to discretise the model.

The finite element results computed are shown in Fig. 9, including the force versus displacement response as well as contours of phase transformation and phase field fracture parameter. Contrarily to the response observed in the previous case study, Section 4.2, a sharp drop in the force versus displacement curve is observed, indicative of brittle fracture with little inelastic dissipation. This can be rationalised by observing the phase transformation contour just before failure, at $u^{\infty}=0.39$ mm; as shown in Fig. 9, the inelastic region is confined to a small area in the close vicinity of the tips of the notches. No inelastic shear bands are observed such that as soon as cracking initiates at the notch tips, the cracks coalescence and unstable cracking is observed. We note that this finding is specific to the boundary value considered, a parametric analysis is needed to characterise the interplay between the different scales at play (notch radius, phase field length scale, sample dimensions) and its implications on fracture stability due to inelastic dissipation.

The capabilities of the modelling framework in capturing fatigue crack growth are demonstrated by simulating cyclic damage in a stent, a paradigmatic application of shape memory materials [69, 70, 71]. This case study also serves to showcase the computational efficiency of the framework and its applicability to the modelling of computationally demanding large-scale 3D boundary value problems.

Stents are small cylindrical tubes that are placed into blood vessels, arteries, or other ducts to hold the structure open. Often, their role is to counteract the effects of vascular diseases that are associated with plaque blockages that hinder fluid flow, see Fig. 10. To deploy the stent, it is usually crimped, placed in a delivery system (e.g., catheter), and finally expanded in-vivo to widen the duct. Nitinol is a popular material choice in stent manufacturing due to its biocompatibility and capacity to expand by recovering its elastic deformation after the constraining delivery system has been removed (superelasticity). However, fatigue resistance is often the limiting design criterion as the stent is subjected to repeated contraction and expansion during the systolic and diastolic cycles. Current fatigue design models for stents commonly use Goodman’s and other empirical methods to estimate the number of cycles to failure by extrapolating from the stress/strain state of the first cycle. We aim here at providing a more mechanistic approach using a phase field fatigue framework that can predict features such as S-N curves or the Paris law as a natural outcome of the model [35].

We assume that the stent has been manufactured using the equiatomic NiTi whose material properties are listed in Table 1. Taking advantage of symmetry along the longitudinal direction, half of the stent is modelled. The stent is subjected to different expansion and compressive pressures by prescribing a radial displacement in the outer surface, $u_{r}$. First, there is a crimp phase where the simulation emulates the compression of the stent inside of the capsule prior to delivery; a radial displacement of $u_{r}=-0.16$ mm is applied. Secondly, the stage of deployment is reproduced by allowing the stent to expand up to $u_{r}=-0.02$ mm, where further expansion is limited by the surrounding duct. Once deployed, the stent is subjected to cyclic loads of amplitude $\Delta u_{r}=-0.04$ mm that simulate the compression and expansion pressures experienced due to the systolic and diastolic cycles. The damage response will be governed by the phase field fatigue model described in Section 2.2.1, with a phase field length scale of $\ell=0.02$ mm and using the monolithic quasi-Newton solution scheme. The magnitude of $\ell$ is at least 2.5 times larger than the characteristic element length. Eight-node brick elements with full integration are used to discretise the geometry, with the model containing more than 7 million degrees-of-freedom (DOFs). Calculations are run in parallel, using 16 cores, with each load increment taking approximately 2 minutes.

The finite element results obtained are shown in Figs. 11, 12 and 13, along with a video that is provided in the online version of this manuscript. Fig. 11 shows the deformed shape of the model during the four stages of the analysis: the crimp and expansion stages of the deployment phase, and the systolic compression and diastolic expansion stages associated with each pressure cycle. Contours of the phase field are shown with the deformed shape, revealing that values of $\phi$ of up to 0.56 are attained during the deployment phase. Thus, a significant amount of damage occurs during deployment but no cracking is observed. This is also the case for the phase transformation; most of it takes place during the deployment phase. The contours of martensitic volume fraction $\xi$ at the end of the crimp-expansion deployment process are shown in Fig. 12.

The $\xi$ contours show that the martensite phase is localised in the edges of the stent struts. This is also the location where cracking initiates. While $\phi$ reaches high values during the deployment stage, 29 cycles of (compression-expansion) systolic-diastolic pressure are needed for cracks to initiate. As shown in Fig. 13, these surface cracks initiate in the regions where phase transformation took place during stent deployment. The cracked region ($\phi=1$) extends with increasing the number of pressure cycles, and after 50 cycles it has extended over a significant part of the stent, including the bridging areas between struts. The evolution of the phase field over 6 cycles is shown in Video 1, provided in the online version of this manuscript. In summary, the example demonstrates that the framework presented here can be used to estimate the lifetime of medical stents for arbitrary geometries, boundary conditions and material properties.