Critical assessment 31: on the modelling of tertiary creep in single-crystal superalloys

ABSTRACT

Modelling of tertiary creep in single-crystal superalloys – operative over the majority of the temperature/stress regime of relevance – is assessed. Traditional, empirical approaches are useful for the prediction of component life; a modern trend is for more physically faithful models which account for microstructure and alloy composition. Calculations are made using different methods, to identify advantages/disadvantages, in an effort to approach best practice. The predictions are tested against experimental data. The creep resistance can be estimated to within approximately one-half of a Larson–Miller parameter, from an input of alloy composition, deformation conditions, and estimates of microstructural parameters, rate laws (and parameters therein), and physical constants. Different merit indices are considered to identify the most appropriate for compositional optimisation.

KEYWORDS:

• Single-crystal nickel-based superalloys
• creep damage
• continuum damage mechanics
• creep modelling
• alloy design

1. Introduction

The nickel-based superalloys – especially when cast in single-crystal form – display outstanding resistance to creep behaviour even at temperatures beyond 1000°C. This justifies their use for critical components needed for jet propulsion and power generation [1–3]. Nevertheless, when one examines typical compositions employed in practice, one finds that only relatively few alloys have entered service. Over time, there has been a trend to improve creep performance via alloying with refractory elements such as tantalum, molybdenum, tungsten, and particularly rhenium, at the expense of reduced corrosion resistance impaired by the lower chromium contents needed to maintain alloy stability, e.g. with respect to the precipitation of topologically close-packed (TCP) phases [4,5].

Why are relatively few alloy compositions available? The main reason is the traditional, overarching emphasis on empirical means of alloy development, which involves many make/test iterations, followed by experimental testing for the purpose of screening, ranking, and subsequent down-selection. This leads to considerable cost but also significant lag in the insertion of new grades into practice [6]. There is no doubt also that the theory and models developed so far have been insufficiently accurate to make a difference to the prevailing situation. Consider the available models for the creep deformation behaviour of this class of alloy, for example, those based upon damage mechanics approaches. The earliest contain many fitting parameters and are rather empirical [7–11]. The latest are more microstructurally explicit and depend upon estimates of microstructural-scale and rate-dependent evolution laws for dislocation content [12–16]. Most can be applied well once a given alloy composition has been fixed, for example, to predict the effects of altering temperature or stress state. So they are useful for component lifing. But it is very rare for such models to be sensitive to changes in material chemistry; it is only with this capability that optimisation of alloy composition can be achieved [17].

The research in the present paper was carried out with the above in mind. We set out to consider the modelling of the creep phenomenon in these materials, with an emphasis on the influence of alloy composition, microstructural architecture, and dislocation evolution. Within this scope, we focus on the tertiary creep deformation mode, operative at intermediate levels of applied stress and temperature, as illustrated in Figure 1. Its main characteristics are that neither precipitate shearing nor rafting is prevalent and that deformation is largely confined to the γ matrix, thus giving rise to a steady increase in creep strain rate with creep strain [1,18,19]. Existing analysis procedures available in the literature are reviewed first. Then, a new merit index which relates alloy composition to creep performance is introduced. By testing against a significant set of experimental data mined from the literature, it is demonstrated that the new approach has advantages over those currently available.

Figure 1. Plot illustrating the regimes of temperature and applied stress under which primary, tertiary, and rafting modes of creep deformation are dominant for the single-crystal superalloy CMSX-4, adapted from reference [1].

2. Background and literature review

Early attempts to describe the evolution of creep damage in engineering materials by da Andrade [20], Dorn [21], Garofalo [22], or Webster [23] were based largely upon empiricism and sought to offer optimal mathematical representations of observed experimental creep curves. Later advances including the 6-θ-projection technique [24] or the Ω-method [25] added more fitting parameters to their formulations. While this facilitated reasonable mathematical descriptions of large creep data sets, only limited further insight into operative creep damage mechanisms was gleaned.

Such shortcomings led to the development of creep continuum damage mechanics (CDM) models [26,27]. These relate the creep strain rate to measurable, external parameters – stress, strain, temperature, and time – and to a number of internal state variables, , which provide a measure for accumulated damage [26]. A constitutive model is then proposed by an equation for the evolution of creep strain rate and by i equations for the evolution of : (1) The origins of creep CDM lie in the attempts of Kachanov [28] and Rabotnov [29] to quantify damage evolution during tertiary creep using a single internal damage parameter. This method was later expanded by Leckie and Hayhurst [30] to also include the effects of multiaxial stress states. Despite these approaches being similar to previous empirical relations, they laid the foundations for physically based CDM, in which damage parameters and damage rates are associated with distinct physical mechanisms and are directly related to microstructural evolution. Dyson [31] identified three broad categories of operating mechanisms based on the kinetics of damage evolution: (1) strain-induced damage, where the damage rate depends on stress, temperature, the level of internal back-stress, H, and the current level of accumulated damage, D; (2) thermally induced damage, depending on temperature and accumulated damage alone; and (3) environmentally induced damage, in which the mean radius of the specimen, R, also plays a role in controlling the damage rate: (2) A comparison of Equations (1(1) ) and (2(2) ) reveals that H and D represent two physically rationalised state variables, and , for which the evolution functions and are derived from theory. Suitable damage parameters and damage rates have been compiled by McLean and Dyson [27] for the most common examples of creep damage: multiplication of mobile dislocations and creep-constrained cavity growth (strain-induced), particle coarsening and depletion of solid solution elements (thermally induced), and internal oxidation or corrosion (environmentally induced). The present work focuses on the first two categories and assumes that environmental effects can be neglected for bulk specimens, as do virtually all other creep CDM models for nickel-based superalloys.

A selection of constitutive equation sets describing tertiary creep in superalloys is presented chronologically in Table 1. A general comparison reveals a persisting lack of agreement on the impact of certain factors on creep resistance. However, the overall approach has remained the same and can be discussed in the context of physically based CDM [34]. For example, an early creep model proposed by Dyson and McLean [7] can be re-written as a CDM formulation involving multiplication of mobile dislocations and particle coarsening. Furthermore, Dyson [31] showed that empirical lifing methods including the Larson–Miller parameter can be readily derived from CDM after making suitable physical assumptions.

Table 1. Models of tertiary creep in superalloys compiled from open literature. Each model is represented by a set of constitutive CDM equations describing creep damage evolution. The column on the right lists parameters determined by fitting models to experimental results and which are not given a clear physical derivation in the respective studies. A general trend towards less fitting parameters over time can be observed.

Study	Year	Constitutive model equations	Fitting parameters
Dyson and McLean [7]	1983
Ion et al. [8]	1986
Dyson and Osgerby [11]	1993
McLean and Dyson [27]	2000
Manonukul et al. [12]	2002
Miodownik et al. [32]	2003
Basoalto et al. [13]	2004
Reed [1]	2006
Dyson [14]	2009
Oruganti [33]	2012
Zhu et al. [15]	2012
Kim et al. [16]	2016

Early creep models [8–11] provide few or no explanations for the values of fitting parameters involved. Although this issue was acknowledged later on [31], most parameters were still obtained by fitting numerical calculations to experiments. The first truly physically based model, relating damage rates directly to alloy composition and microstructure, was that of Manonukul et al. [12]. It was followed by other similar formulations with certain adjustments [13–15,35]. For example, Coakley et al. [36,37] expanded a previous model by Dyson [14] to account for multimodal particle distributions. A notable recent entry is the work of Kim et al. [16], who followed a similar approach to the present study and used thermodynamic calculations to inform their creep model with regard to microstructure. It should be pointed out that the list in Table 1 is only exemplary, and that many other models have been proposed, some of which include mechanistic expressions too complex to be reproduced above in short [35,38–45].

When comparing physically based models for intermediate-temperature creep, several common traits within the expressions for creep strain rate become clear. First, creep resistance increases with higher volume fraction, , and lower effective diffusivity, . Second, the internal back-stress, usually depending on particle spacing, , has a profound significance on creep rate evolution. Third, particle coarsening plays an important role in affecting particle spacing, and in turn internal stresses. Fourth, the choice of values defining initial microstructure, including initial dislocation density, , and initial precipitate radius, , influence model predictions significantly. Finally, the same is also the case for functions derived to describe parameters like effective diffusivity, volume fraction, or an invoked constant damage multiplication parameter, or , used to quantify the increase in mobile dislocation density. The latter two points have been clearly illustrated in the work of Zhu et al. [15], who observed significant changes in simulated creep curves with varying values of volume fraction, particle size, initial dislocation density, and dislocation multiplication parameters. It can be expected that other models will behave similarly. This is problematic in the context of developing a model suitable for a wide range of alloys, as constant values rationalised for one particular case may not apply to others.

A critique on the usefulness of the different models is challenging for a number of reasons. Oruganti [33,46] pointed out that cases of poor correspondence between experiments and model predictions are often attributed to gaps in understanding damage mechanisms rather than to the basic constitutive creep laws used. As such, the critical assessment of these laws is of importance. For example, efforts to describe creep over a wide range of conditions led in the 1940s to the use of exponential or hyperbolic sine laws [47]. These were later rationalised theoretically by Friedel [48] and Hirth and Lothe [49] in their seminal works on dislocation theory. An advantage of expressions is their ability to capture the strong influence of applied stress on creep life. However, when used in some models above, they give rise to an overestimation of the role of particle spacing, which is not backed up by experiments [27,33]. Furthermore, the use of constant stress exponents within power-law descriptions has been recently criticised as too simplistic and not supported by microstructural and physical observations [50]. The same can also be stated in the case of internal back-stresses, often associated with certain strengthening mechanisms and assumed to evolve over time [51,52]. Attempts at addressing such issues with the introduction of effective – as opposed to directly measured – creep parameters do not shed more light on intrinsic damage accumulation [53,54]. Lastly, damage multiplication parameters, which relate the increase in mobile dislocations to the respective strain rate, have been assumed to depend linearly on initial dislocation density and to remain constant [13,14]. No clear explanation for this choice exists, and recent discrete dislocation dynamics (DDD) modelling efforts revealed that severe fluctuations in the value of can be expected during creep [55]. Such shortcomings provided an incentive for the development of a new modelling approach, complex enough to capture key influences of chemistry and microstructure on alloy properties, but simple enough to be applied in alloy design studies.

3. Modelling of tertiary creep phenomenon

3.1. General modelling approach

The main procedural steps and computational aspects used in this paper build upon those of the Alloys-By-Design (ABD) method [6] and are illustrated in the process flow diagram of Figure 2. Necessary inputs include alloy chemistry and creep conditions – applied stress along the technically relevant crystal direction and temperature. These inputs can be precisely specified for known commercial alloys or standard testing conditions. Alternatively, a large design space can be generated by setting ranges and increments for each chemical element and can then be investigated to isolate compositions predicted to possess optimal property trade-offs, following references [56,57].

Figure 2. Process flow diagram representing the main steps taken to model precipitate coarsening and creep strain accumulation.

Alloy composition and temperature serve as inputs for CALPHAD thermodynamic modelling which estimates key microstructural parameters. In the present work, version 2016b of the Thermo-Calc software package was used along with the TCNi8 thermodynamic database and the MOBNi3 mobility database. Computed values include the interfacial energy between γ and , , the partitioning coefficients of different elements, and phase volume fractions. Of particular interest are the fraction, , and – for less stable alloys – the fractions of deleterious TCP phases. In combination with further models and data from the literature, Thermo-Calc predictions are used to estimate other important parameters such as lattice misfit, δ, and the coarsening rate coefficient, k. As noted by the use of dashed arrows, the current approach can be extended with other simulation packages. This includes using TC-PRISMA to analyse effects of chemistry and heat treatments on nucleation and growth of precipitates, or using DICTRA to estimate effective diffusivity parameters.

In the final steps, sets of constitutive model equations are integrated for each alloy from time zero using standard finite difference methods. The results are predictions regarding (1) the coarsening of precipitates during high-temperature exposure and (2) the evolution of creep strain. Further user-defined merit indices – such as density-corrected creep strength – can be applied within ABD studies.

A novel addition is incorporating data from correlative high-resolution EBSD (HR-EBSD) and electron channelling contrast imaging under controlled diffraction conditions (cECCI) studies on specimens tested until pre-defined levels of creep strain. These provide both a quantitative and a qualitative evaluation of dislocation structure evolution during creep and lead to a more experimentally grounded description of damage through dislocation accumulation.

3.2. Modelling of superalloy microstructure and physical properties

Most creep models in the literature focus on reproducing experimental results for a few commercial alloys, and so they apply known values of precipitate size, elastic modulus, or lattice misfit. However, for alloy design purposes, changes of these properties with chemistry must be considered in further detail. In certain cases, this cannot be achieved accurately with current CALPHAD methods, and pragmatic estimates are needed.

3.2.1. particle distribution

Essential features of the microstructure being modelled are shown in Figure 3. For modern single-crystal superalloys, it consists of a high volume fraction of the ordered precipitate phase residing within the face-centred cubic γ matrix. Precipitates are modelled as cubes with side length d and interparticle spacing w. A simple, unimodal distribution is shown below; however, the model accounts for variations in size by assuming a normal size distribution for d, with a mean value of 450 nm and a standard deviation of 150 nm. As the constitutive creep equations require discrete values, this distribution is approximated by a histogram containing 5 bins with an equal bin width of 100 nm.

Figure 3. Assumed evolution of damage during tertiary creep in single-crystal superalloys. The main operating mechanisms are (1) an increase in mobile dislocation density in γ matrix channels, ρ, and (2) coarsening of precipitates, i.e. an increase of the average cube side length, d, and of the interparticle spacing, w.

The mean value of 450 nm has been noted to lead to optimal creep properties in typical modern superalloys with a low, negative lattice misfit at high temperatures [58,59]. It is assumed here that it can be attained for any particular composition through appropriate heat treatments. However, it would be insufficient to only consider mean particle size – as done in most creep models – since it has been shown that a homogeneous microstructure, with evenly distributed particles oriented along the crystal direction, significantly improves creep strength [60]. The standard deviation of 150 nm was chosen based on particle size measurements carried out on recently developed ABD alloys [56] and is also assumed to remain constant for all investigated compositions. For regions with constant volume fraction and size, it follows from geometry that , and that .

3.2.2. Elastic properties

Influences of chemistry on elastic properties have been commonly estimated separately for each phase using ‘rule-of-mixtures’ approaches, with the moduli of individual elements being scaled for temperature [16,61]. For example, Kim et al. [16] calculated the temperature dependence of shear modulus, G, as: (3) Here, factors represent changes in the shear moduli of pure nickel and of pure Ni₃Al caused by an amount of alloying element i. For the phase, different changes in modulus depending on the location of alloying elements on Ni or on Al sites are taken into account. Partitioning coefficients and site fractions can be obtained from CALPHAD calculations. Ultimately, Kim et al. derived optimal values for , , and for a range of elements by calibrating Equation (3(3) ) to match experimental data.

While this approach shows promise for general alloy design purposes, it assumes isotropic elastic properties and can only be applied to polycrystalline materials. In contrast, early studies of elastic properties in directionally solidified superalloys showed a strong dependence on crystal direction [62,63]. As discussed in more detail by Reed [1], this stems from the anisotropic elasticity of pure nickel in single-crystal form, and must be considered when modelling mechanical properties of new single-crystal superalloys.

Siebörger et al. [64] carried out extensive measurements of the temperature dependence of Young's modulus, shear modulus, and Poisson's ratio of superalloy CMSX-4 and its isolated γ and phases. Based on tests in near and directions, these authors derived linear temperature dependences for , , and between room temperature and 1000°C. These data can be used to calculate stiffness matrix components , , and at a certain temperature. Lastly, elastic properties along the active slip planes during loading along the tensile axis can be obtained.

Morrell et al. [65] compared the relative effects of temperature, crystal orientation, and chemical composition on the elastic moduli of several commercial single-crystal superalloys. These authors concluded that the effect of composition on stiffness tensor values is a minor one. This is in agreement with the observations of Dye et al. [66], who compiled experimental data on stiffness matrix components for a range of superalloys and compared these results to the levels of solid solution strengtheners (Cr+Mo+Re+W) and of -formers (Al+Ta+Ti) in each alloy. No clear trends could be distinguished for solid solution strengtheners, while titanium in particular was noted to increase the stiffness of . Considering the results of Morrell et al. and Dye et al. pragmatically, it can be expected that deviations from the ideal orientation during casting will have a much more significant effect than chemistry. As such, in the present work, the temperature-dependent results obtained by Siebörger et al. for CMSX-4 are applied generally to all compositions.

3.2.3. γ/ lattice misfit

Interfaces between matrix and precipitates are assumed to be coherent in the present study, and misfit stresses are not taken into account. As lattice misfit is often cited as a key parameter for alloy design and improved creep strength, reasons for this omission are discussed in the following.

First, there is no clear consensus regarding either the optimal magnitude, or the optimal sign of lattice misfit for increased creep resistance [67]. On the one hand, it seems beneficial to keep the value close to zero, in order to obtain a low coarsening rate coefficient and homogeneous, cuboidal precipitates [1,68–70]. On the other hand, it has been suggested that a high, usually negative, mismatch will promote dense dislocation networks at γ/ interfaces, ultimately leading to increased protection against precipitate shearing [71–79].

Second, regression analysis studies carried out by Harada et al. [80–83] revealed a decreasing influence of misfit and an increasing influence of solid solution strengthening on creep properties as testing temperature is reduced. While dense interfacial dislocation networks appear important in the rafting regime, they do not play a main role for tertiary creep.

Third, and of most importance from a practical perspective, influences of chemistry and temperature on lattice misfit cannot be accurately predicted by current methods, as seen in Figure 4. With respect to misfit calculations for alloy design, most authors including Caron [84], Reed [1,6], and Rettig [70] follow the approach proposed by Watanabe and Kuno [85] to calculate lattice misfit, δ, and lattice constants for γ and at room temperature as: (4) where is the lattice constant of pure nickel or pure Ni₃Al, and and are the Vegard coefficient and the atomic fraction of chemical element i in phase j. Following Caron [84], and are taken as 3.524 Å and 3.57 Å, respectively. Values for the Vegard coefficients of different elements can be found in the works of Mishima et al. [86], Nash [87], Kablov et al. [88], or Reed [1]. Phase compositions can be estimated using CALPHAD methods. While reported values for Vegard coefficients are in reasonable agreement, there are significant differences between studies for certain elements, leading to discrepancies between misfit results based on different datasets.

Figure 4. Comparison of changes in lattice misfit with temperature from (a) experimental XRD measurements [92], (b) empirical calculations of lattice constants [70], (c) empirical calculations of lattice constants and thermal expansions coefficients [90], and (d) CALPHAD calculations of molar volumes.

The influence of temperature on lattice constants for pure nickel and Ni₃Al has been reported in the work of Kamara et al. [89] and Reed et al. [6], and is also taken into consideration in the alloy design studies cited above. However, chemical additions will influence not only lattice constants at room temperature, but also thermal expansion coefficients. The most comprehensive studies covering this topic are the empirical model of Sung and Poirier [90], based on experimental data for a range of superalloys, and the phase mixture model of Karunaratne et al. [91]. Both approaches show good agreement with selected literature data, however, prove difficult to be applied more broadly.

A promising and straightforward approach is to obtain lattice constants for the γ and phases from molar volumes calculated by CALPHAD, which accounts for both chemical composition and temperature. As these phases contain four atoms in each unit cell, lattice constants can be obtained as , where is the calculated molar volume of phase j and is the Avogadro constant.

Unfortunately, when compared to the experimental work of Pyczak et al. [92], none of the approaches discussed above can accurately predict changes in lattice misfit with temperature, as shown in Figure 4. Pyczak et al. carried out X-ray diffraction (XRD) measurements on a range of rhenium-containing and rhenium-free superalloys and noted an increasingly negative mismatch at higher temperatures for all compositions (Figure 4(a)). This qualitative trend is reproduced correctly by the misfit formula used by Rettig et al. [70] for alloy design, although quantitative values differ markedly (Figure 4(b)). Misfit values obtained based on the empirical model of Sung and Poirier [90] (Figure 4(c)) or based on CALPHAD calculations of molar volumes (Figure 4(d)) deviate even more from experiments. Lastly, it is interesting to note that all three modelling approaches in Figure 4(b)–(d) largely predict the same ‘descending order’ of lattice misfit from IN792DS to ExAl8, although this is not in good agreement with XRD results.

3.3. Modelling of damage evolution during tertiary creep

3.3.1. Kinetic creep equation

The derivation of the kinetic creep equation follows the considerations of Dyson [14] and Zhu et al. [15], who identify dislocation accumulation in γ matrix channels and coarsening of precipitates as the main sources of damage. It is acknowledged that strain-induced nucleation and growth of cavities and microcracks will occur in the later stages of creep [18,19]; however, these additional damage mechanisms are omitted in the present model, which focuses on the initial deformation regime up to several percent plastic strain.

In the undeformed state, only a few dislocations are located within γ channels, as represented schematically in Figure 3. Recent cECCI work by Ram et al. [93] showed conclusively that such isolated dislocations are the main dislocation sources during tertiary creep. Under an applied uniaxial stress, some of them will become active via Frank-Read [94] or double cross-slip mechanisms [95] once a critical shear stress has been reached locally. Low-angle grain boundaries (LAGB) formed as crystal defects during casting [96] also give rise to dislocations; however, these were observed to only percolate a few micrometers away, meaning that LAGBs are not major sources as previously assumed [93].

Newly generated dislocations glide until reaching the nearest particle, where they become arrested at γ/ interfaces. Depending on the local geometry and on the sign of uniaxial stress, a series of climb events through vacancy absorption or vacancy emission are required to overcome precipitates [14]. For the usual case of an applied tensile load, vacancy absorption will take place in ‘vertical’ channels and emission in ‘horizontal’ ones, leading to creep deformation progressing through combined climb/glide events. This situation has been theoretically assessed by Friedel [48] and Hirth and Lothe [49], who derived expressions for the jump frequencies of vacancy-absorbing and vacancy-emitting jogs. These can in turn be combined to yield a net escape frequency under the premise of an equal number of absorption and emission events. The particular rate at which trapped dislocations are released constitutes the core of previous creep models by Dyson [14] and Zhu et al. [15]. However, a linear dependency between creep strain rate and is debatable, as discussed in Section 2.

In the present work, a different approach is taken by considering relations between (1) effective forces acting on dislocations under the influence of solute drag and (2) average dislocation velocities, for both the climb and the glide step. This approach builds on insights gained from recent DDD studies [55,97,98] and also considers novel experimental work by ur Rehman et al. [61] on quantifying climb and glide drag coefficients in the γ matrix. It is acknowledged that not all theoretical contributions towards the forces acting on dislocation segments are considered [98]. Some sources such as the line tension which minimises the dislocation core energy, the elastic force due to the stress field of the dislocation network, or the osmotic force due to vacancy depletion or supersaturation are too complex in nature to be included in the present model. Furthermore, following the simplified approach proposed by Takeuchi and Argon [99], the same ‘effective stress’ acting on dislocations is assumed for both the climb and the glide step.

The average dislocation velocity can be expressed as , where τ is the effective shear stress, b is the length of the Burgers vector, and B is the drag or friction coefficient [100]. It can be further related to the microscopic shear strain rate, , via the Orowan equation, , where is the total density of mobile dislocations. A few points must be taken into consideration. First, the total density is composed of the density of climbing dislocations at particle/matrix interfaces, , and the density of gliding dislocations within matrix channels, . Second, a distinction must also be made between the average velocities of climbing and gliding dislocations, and , owing to different drag coefficients, and . Third, the standard Orowan equation cannot be readily applied to precipitate-strengthened materials like superalloys, since to a large extent, deformation is localised within the matrix and precipitates only deform elastically during tertiary creep. This point is addressed by Dyson [14] by considering the necessary load transfer from matrix to particles as plastic creep strain accumulates. Assuming that the elastic moduli of matrix and precipitates are equal, a simple relation between the measured global shear strain rate, , and the microscopic shear strain rate within the matrix, , is given by . The premise that no precipitate shearing occurs is not entirely true; however, TEM and ECCI studies confirm that it is not significant during tertiary creep [18,93]. Finally, following the considerations of Ghosh et al. [101], eight slip systems will be active during symmetrical loading along the tensile axis and will have equal Schmid factors of . Thus, the macroscopically measured creep strain rate, , can be expressed as . Overall, we obtain: (5) Attention should now be given to the effective operating shear stress, τ, resulting from externally applied stress on the one hand, and strengthening mechanisms on the other. Following the considerations of Pollock and Field [18], creep strain rate in the tertiary regime is obtained as: (6) where is a material constant, is the applied stress, and , , and are the respective contributions from Orowan strengthening, solid solution strengthening, and dislocation resistance. Some models also include an additional stress component arising from lattice misfit, , which is added to the applied stress. As most commercial alloys exhibit negative misfit at high temperatures, this superposition leads to significantly increased deformation in ‘horizontal’ γ matrix channels, perpendicular to the loading direction, and only limited creep in ‘vertical’ channels, parallel to loading [102]. Considering the aforementioned issues regarding ABD calculations of misfit, it is both difficult and open to doubt to incorporate misfit stresses in a model used for alloy design, as incorrect predictions of misfit – and therefore of misfit-induced stress – could lead to faulty conclusions and design rules.

In Equation (6(6) ), the main source of strengthening is the local Orowan resistance of matrix channels, which must be overcome in order to bow dislocations between precipitates. The resolved shear stress required can be expressed as: (7) where G is the shear modulus, l is the distance between precipitates in direction, and w is the γ channel width. A small contribution is also provided by dislocation passing stresses and can be approximated for the case of random interactions between dislocations as: (8) A value of 0.1 is chosen for α, in accordance with references [18,103]. A creep stress exponent, n, is not included in the present model, as such exponents have been shown to only remain constant over a restricted stress and temperature range [31]. As a correction, Dyson [31] derived a more advanced expression for n, in relation to empirical hyperbolic sine functions, as , with being a function of temperature, volume fraction, and Orowan resistance stress. Calculated values for n would normally range between 1 and 2 for nickel-based superalloys, in agreement with the expected dominating mechanism of diffusional climb. Nonetheless, such an expression ultimately does little towards improving the basic understanding of creep deformation [50].

The constraint factor in Equation (6(6) ) accounts for the additional frictional resistance imposed by non-deforming particles on the matrix, which results in large pressure gradients in γ channels. This factor is discussed in detail by Pollock and Argon [103], who offer an analytical approximation as a function of precipitate side length, d, and channel width, w. Assuming an ideal, unimodal particle distribution like in Figure 3, the constraint factor can be expressed as depending only on the volume fraction: (9) This expression is incorporated readily into the models of Basoalto et al. [13] and Zhu et al. [15]. However, Pollock and Argon [103] debated its general validity, as it leads to values of tensile creep resistance far in excess of experimental estimates. Instead, it appears that as γ channels become populated by a knitted network of dislocations, many of which are lying along γ/ interfaces, some additional sliding will become possible [103]. As such, Pollock and Field put forward a value between 1 and 2 for ζ, depending on the size and shape of precipitates [18]. Van Sluytman and Pollock [104] discussed how such a factor might vary with morphology, and in turn with lattice misfit. They predict optimal values can be reached in the range of . However, since accurate misfit predictions are again required, it is difficult to incorporate a shape- and misfit-dependent factor into the present model.

Furthermore, this resistance from the build-up of pressure gradients due to frictional drag will evolve during the creep test. Dyson [14,31] took ζ to increase from unity to a maximum value of due to stress redistribution around particles and proposed a stochastic approach to describe the evolution of internal stress, , as: (10) where is the upper bound collapse value. An analogue relation to Equation (10(10) ) is constructed to describe the evolution of ζ with creep strain rate, and the upper bound value of is chosen as 1.75. This is slightly lower than the initial, ideal estimate of Pollock and Argon from Equation (9(9) ) and the simple expression of proposed by Dyson [31], in order to account for imperfect particle shapes and reduced resistance to sliding at particle interfaces.

Equation (5(5) ) for creep strain rate can now be expressed as: (11)

3.3.2. Effective climb and glide drag coefficients

As a consequence of the substantial levels of solid solution strengthening in modern superalloys – normally maximised up to the limit of TCP phase precipitation [70] – dislocations experience significant drag forces and move through γ channels in a viscous manner [103]. Therefore, the values of the effective climb and glide drag coefficients in Equation (11(11) ), and , play a key role in determining creep strain rate.

Early estimates come from the theoretical work of Cottrell and Jawson [105] and Hirth and Lothe [49], who considered gliding dislocations surrounded by an atmosphere of solute atoms at high temperatures. Takeuchi and Argon [99,106] expanded on these results and derived expressions for modelling the effects of viscous drag on climbing and on gliding edge dislocations. In more recent studies, ur Rehman et al. [61] applied the model of Takeuchi and Argon to assess climb and glide resistances in single-crystal binary Ni–X (X=Re,Ta,W) solid solution alloys, while Hafez Haghighat et al. [55] and Liu et al. [98] used similar descriptions as part of DDD simulations of tertiary creep in single-crystal superalloys. The authors of the latter study expressed the drag coefficient for climb as: (12) where is the Boltzmann constant, is the diffusion coefficient of the solute, is the solute concentration at a theoretically infinite distance from the dislocation, and Ω is the atomic volume of the solvent. Due to the inverse relation between and , additions of slow-diffusing elements – particularly molybdenum, rhenium, tantalum, and tungsten for superalloys – will have a profound effect on generating increased drag. Similarly, the drag coefficient for glide was evaluated as [61]: (13) with factors as defined above. While solute diffusivity again plays an important role, it is supplemented by misfit effects due to differences in the atomic radii of solute and solvent. These are expressed by the size interaction parameter, β, which can be obtained as [99]: (14) The term , present in both expressions, is used to define the interaction between dislocations and the solute environment, with and being the inner and the outer cut-off radii. Following the considerations of Hirth and Lothe [49], is equal to the length of the Burgers vector, b, while depends on the velocity and the diffusivity of the solute as . However, when this ratio becomes larger than the average distance between dislocations in the crystal – which is the case for densely packed dislocations in matrix channels – it no longer represents a meaningful estimate and must be replaced by an appropriate distance, R [49]. Pollock and Argon [103] used the value of γ channel width, w, as an estimate for R. Considering microstructural coarsening during creep, the term will increase gradually towards a constant value. As a simplification, we assume this value to be effective throughout the duration of creep tests. Taking w equal to 70 nm as a representative value for the coarsened state [107] yields , similar to the constant factor of 4.3 first derived by Cottrell and Jawson [105].

So far, only the influence of an individual solute species has been considered. In order to estimate total, effective drag coefficients, we follow the considerations of Reed et al. [6] with regard to dislocation pinning by different types of solute atoms. These authors defined an effective time to circumvent a unit length of pinning atoms as: (15) Terms denote the site densities of pinning solute atoms met by the leading segment of the dislocation, while are the characteristic times of escape, depending on the relative pinning strength of the solutes. The time of escape will scale inversely with average dislocation velocity, and therefore linearly with the drag coefficient, as . As such, total, effective drag coefficients for climb and for glide can be estimated as: (16) and (17) The factor , with and being the mole fraction and the interdiffusion coefficient of solute i, respectively, is identical to the merit index for creep resistance proposed by Reed et al. [6]. For single-crystal superalloys with a similar microstructural architecture – in terms of lattice misfit and volume fraction – a clear correlation between the Larson–Miller parameter and this merit index has been established [6,56,57]. This correlation was strengthened when the calculated solute mole fraction in the γ matrix, , was used instead of the mean composition. The same is done in the present study, as it is reasonable to calculate drag coefficients based on the composition of the γ phase.

Interdiffusion coefficients, , are calculated as . Pre-exponential factors, , and activation energies pertaining to the interdiffusion process, , have normally been obtained experimentally using Ni–X diffusion couples. The ample work of Karunaratne et al. [108–111] covered all relevant chemical additions in superalloys with the exception of chromium and cobalt, which had been assessed previously [112,113]. A critical review reveals significant discrepancies between experimental results in different reports and makes it difficult to choose ‘optimal’ values for estimating effective diffusivity [113]. For alloy design purposes, a CALPHAD-based approach using a multicomponent diffusion mobility database, as proposed by Campbell et al. [114], appears more promising. A comparison between experimental and CALPHAD-based results, along with their respective sources, are shown in Table 2. While most values are in good agreement, it is clear that using experimental rather than CALPHAD-based data will give more weight to tungsten and rhenium additions. It can be argued that CALPHAD-based data will help minimise experimental scatter by considering several sources, as such, values in the ‘CALPHAD’ column are used in the present study.

Table 2. Pre-exponential factors and activation energies for solute diffusion in nickel based on diffusion couple experiments and CALPHAD predictions.

	Experiment			CALPHAD
Element	[m^2 s^−1]	[kJ mol^−1]	Reference	[m^2 s^−1]	[kJ mol^−1]	Reference
Cr	5.20×10^−4	289	[113]	4.33×10^−4	287	[114]
Co	1.80×10^−4	282	[113]	2.94×10^−4	284.2	[114]
Mo	1.15×10^−4	281.3	[111]	1.19×10^−4	281.8	[115]
Re	8.20×10^−7	255	[108]	3.28×10^−6	278.8	[114]
Ru	2.48×10^−4	304.4	[110]	2.52×10^−4	304.5	[116]
W	8.00×10^−6	264	[108]	2.79×10^−5	282.1	[114]
Al	2.30×10^−4	265	[109]	7.52×10^−4	284	[114]
Ti	2.30×10^−4	269	[109]	9.06×10^−5	256.9	[114]
Ta	2.19×10^−5	251	[108]	6.71×10^−5	267.7	[114]

Finally, it is instructive to estimate the ratio at operating temperatures for successive generations of superalloys. This serves both as a sanity check, by comparing these results to usual assumptions in DDD studies, and as a critical assessment of the effects of solid solution strengthening in the present model. Following Zhu et al. [117], the alloys SRR99, CMSX-4, TMS-75, and EPM-102 are taken as representative of first- to fourth-generation superalloys. Their nominal chemical compositions are shown in Table 3.

Table 3. Nominal chemical compositions of exemplary first- to fourth-generation single-crystal superalloys in wt.%.

Generation	Alloy	Cr	Co	Mo	Re	Ru	W	Al	Ti	Ta	Hf	Ni	Reference
First	SRR99	8.0	5.0	–	–	–	10.0	5.5	2.2	3.0	–	Bal	[118]
Second	CMSX-4	6.5	9.0	0.6	3.0	–	6.0	5.6	1.0	6.5	0.1	Bal	[119]
Third	TMS-75	3.0	12.0	2.0	5.0	–	6.0	6.0	–	6.0	0.1	Bal	[79]
Fourth	EPM-102	2.0	16.5	2.0	6.0	3.0	6.0	5.55	–	8.25	0.15	Bal	[120]

Results in Figure 5(a) show an increase in resistance to dislocation motion in higher-generation superalloys. This trend is related to the rising amounts of slow-diffusing elements, particularly molybdenum, tungsten, and rhenium. However, it is important to consider that although these elements partition preferentially to the γ phase, actual values of the partitioning coefficients , and consequently effective levels of solution strengthening, may differ. Since the climbing step controls the rate of creep strain accumulation, it is of particular importance to increase rather than . The plot in Figure 5(b) reveals to what extent this is achieved for different superalloys. CMSX-4 provides the highest relative level of solution strengthening, benefiting from a favourable partitioning of rhenium to γ. In the case of EPM-102, ruthenium has been added specifically to increase and limit TCP phase precipitation [120]. This obvious trade-off brings to mind the difficulty of developing new superalloys with improved specific creep strength, as discussed in the introduction. The range of observed ratios is in good agreement with the commonly assumed constant value of 100 from recent DDD studies [55,61,98].

Figure 5. (a) Estimated drag coefficients for dislocation climb and glide in matrix channels. (b) Estimated ratio of .

3.3.3. Mobile dislocation density

The first type of damage accumulation considered within Equation (11(11) ) is the strain-induced multiplication of mobile dislocations. As discussed in Section 2, most phenomenological models for tertiary creep assume a linear relation for the evolution of mobile dislocation density with increasing creep strain: (18) where the constant dislocation multiplication parameter, , is expressed as a multiple of initial dislocation density, , as . Depending on the study, the assumed value for K is 50 [9], 300 [14], 67 [13], or 100 [15,31]. A first matter of contention is that such differences in the value of K, albeit relatively small, result in significant changes of predicted time to 1% creep strain [15]. Furthermore, recent DDD studies of tertiary creep deformation in superalloys predict a non-constant value for [55,93], as illustrated in Figure 6. In this graph, DDD simulations are compared to HR-EBSD and cECCI measurements for the second-generation, rhenium-containing superalloy ABD-1 and the first-generation, rhenium-free superalloy ABD-3 under deformation conditions of 900°C and 450 MPa. Results consistently show a strong initial increase in at low plastic strain values, followed by an exponential decrease as more creep strain accumulates [55,93]. Furthermore, despite the differences in volume fraction and level of solid solution strengthening, similar trends are observed for both alloys. This change in dislocation density with creep strain can be described using the fundamental equation proposed by Gilman [121] for the dependence of mobile dislocation density on the plastic shear strain in crystalline solids: (19) Here, M is a multiplication parameter leading to an initial rise in mobile dislocation density, and φ is an attrition parameter which expresses the exponential decay over time. Equation (19(19) ) was used successfully to describe the creep response of single-crystal CMSX-4 specimens in the tertiary and rafting regimes [122,123]. Based on values obtained from fitting to experimental results, these authors concluded that M increases exponentially with applied stress and that φ depends on crystal orientation and appears to scale inversely with applied stress. Approaches for estimating M and φ are discussed in more detail in Sections 4.1 and 4.2.

Figure 6. Evolution of dislocation density with plastic creep strain for the single-crystal superalloys ABD-1 and ABD-3 deformed at 900°C and 450 MPa. DDD simulation results are compared to experimental values from HR-EBSD and cECCI analyses and show good overall agreement.

For small strain values – as is the case for common design targets of nickel-based superalloys – it follows that: (20)

As the present model differentiates between climbing and gliding dislocations, their mobile densities must be treated separately. Following the considerations of Dyson [14], slightly modified by Zhu et al. [15] for cuboidal rather than spherical precipitates, the change in the net rate of climbing dislocations can be expressed as: (21) with all terms as defined previously. This change is obtained as the sum of the rate at which trapped climbing dislocations are released from precipitate interfaces and the rate at which gliding dislocations become temporarily arrested at the next particle. The half of the cuboid side length, d/2, and the channel width, w, are assumed to be the average distances that dislocations must overcome in each step [15]. Furthermore, the likelihood of climbing dislocations escaping depends on their probability of encountering a particle in the first place, which is taken to be equal to [14]. As climbing and gliding dislocation densities are codependent, the change in the net rate of gliding dislocations can be obtained as: (22) Manonukul et al. [12] discussed the inclusion of two further terms in Equation (22(22) ) to account for changes in the net rate of pinned dislocations, , arising from individual dislocations being immobilised by the dislocation network and from thermally activated network recovery. However, as argued by these authors, the rate of climbing dislocations can be assumed to dominate the rate of pinning, , therefore can be neglected.

Finally, initial values must be chosen for and . Considering ECCI and TEM results, it can be assumed that all mobile dislocations will begin by gliding through γ channels once an external load has been applied. As such, we take and . An often cited value for in the literature has been 10¹⁰ m⁻² [12–15], though none of these authors provide a clear source for this estimate. Similarly, Preußner et al. [124] chose an initial value of 5×10¹⁰ m⁻². An early TEM study of dislocation densities in cast superalloys does indeed mention densities as low as 10¹⁰ m⁻² for the undeformed state [125], and this value is also in accordance with the estimate of Hull and Bacon [100] for as-cast crystals. Nonetheless, as the choice of this value will severely impact model results [15], and as the material undergoes further processing and machining steps prior to testing, a careful consideration is essential. Recent HR-EBSD and cECCI studies [93] yielded a dislocation density of 5×10¹¹ m⁻² in the undeformed regions of machined creep specimens. As such, in the present work, a figure of 10¹¹ m⁻² is considered more representative of the initial microstructure prior to creep testing.

3.3.4. precipitate coarsening

The second type of damage considered is the thermally induced coarsening of precipitates. The fundamental equation for competitive particle growth can be expressed as: (23) Here, and represent average particle sizes at times t and 0. The exponent n is characteristic for the dominating coarsening mechanism and takes a value of 2 for interface-controlled reactions and of 3 for diffusion-controlled reactions. The latter case corresponds to the widely applied Lifshitz–Slyozow–Wagner (LSW) theory [126,127] for diffusion-controlled Ostwald ripening. Coarsening in modern nickel-based superalloys can be accurately described by Equation (23(23) ), with n=3 and , where d is the average side length of cuboids [128]. The key parameter defining coarsening rate is k, and the way in which it is estimated remains a matter of contention between different theories. An ample discussion of this point can be found in the reviews by Baldan [128,129] on the progress in Ostwald ripening theories and on their application to precipitates in nickel-based superalloys. A similar analysis is beyond the scope of the present work, nonetheless, it is important to carefully choose the method of calculating k, as particle coarsening and the corresponding increase in particle spacing will have a significant effect on creep resistance.

All coarsening theories developed for nickel-based superalloys express k as having the form: (24) Here, is a dimensionless proportionality factor, is the effective diffusion coefficient, is the matrix concentration of the constitutive particle element, is the interfacial energy between matrix and particles, and Ω is the molar volume of the precipitate. Different methods of estimating , , , and Ω have been proposed in the literature.

The LSW theory assumes a negligible precipitate volume fraction and obtains from theoretical considerations . Later studies argued that k should scale with and provided different analytical or numerical approaches to estimate . McLean [130] compared the LSW theory with the modified LSW (MLSW) theory of Ardell et al. [131–133], the Lifshitz–Slyozov encounter modified (LSEM) theory of Davies et al. [134], and the Brailsford–Wynblatt (BW) theory [135] and compiled a table of values for different volume fractions, reproduced in Table 4. Six additional theories, each yielding somewhat different results, have been reviewed and discussed by Baldan [129]. Considering all these discrepancies, it can be concluded that no single theory is optimal and that each will be more or less effective depending on the alloy system and the coarsening conditions involved. Furthermore, while numerical approaches provide certain improvements, for example, by reproducing realistic particle size distributions, they also require extensive computational resources, which remains prohibitive for alloy design purposes.

Table 4. Comparison of proportionality factors estimated by different coarsening theories.

(%)	LSW [126,127]	MLSW [132]	LSEM [134]	BW [135]
2	0.89	2.37	1.02	1.16
5	0.89	3.49	1.09	1.44
10	0.89	5.06	1.17	1.75
20	0.89	7.90	1.30	2.19
40	0.89	13.25	1.54	2.98
60	0.89	18.44	1.84	3.69

In the present work, four analytical methods of calculating the coarsening rate coefficient k are compared to experimental results. These are the standard LSW theory, the MLSW theory, the LSW-type theory of Li et al. [136], and the MLSW-type theory of Ai et al. [137]. The latter two methods have been developed primarily to help guide the development and design of new nickel-based superalloys and include estimates of , , , , and Ω depending on alloy composition and temperature. For the LSW and the MLSW theories, these factors must still be evaluated. Following Ai et al. [137], for the MLSW theory is estimated with a polynomial as . For both theories, is obtained following Reed et al. [6] as , analogue to the calculations described in the section on effective drag coefficients. represents the matrix concentration of the rate-controlling element. According to McLean [130], aluminium is likely to be rate-controlling for both Ni₃Al and Ni₃(Al,Ti). The increasing amounts of slower-diffusing tantalum in modern alloys may, however, affect for Ni₃(Al,Ti,Ta). As such, is taken as the sum of the matrix concentrations of aluminium, titanium, and tantalum, following Li et al. [136]. Lastly, and Ω are estimated from CALPHAD calculations using Thermo-Calc and the associated databases TCNi8 and MOBNi3.

Experimental data for testing the models were selected to cover the first three superalloy generations and include results on first-generation CMSX-2 exposed between 800 and 1050°C [138], on second-generation CMSX-4 exposed between 850 and 1000°C [107], and on third-generation CMSX-10 exposed at 950 and 1050°C [139]. This data set allows to check whether theories correctly reproduce the beneficial effects of increased rhenium additions, as this slow-diffusing element has been shown to play a key role in reducing , and in turn k [140,141]. Results in Figure 7 show that the MLSW theory and the approaches of Ai et al. and Li et al. increasingly overestimate k. This observation is in accordance with the critical review by McLean [130], who concluded that neither particle shape, nor volume fraction appear to have an obvious effect on k, implying that the factor can be largely neglected. The standard LSW theory slightly underestimates k. As such, a slight correction as provided by the LSEM theory would be beneficial but remains challenging to implement due to the associated computational cost.

Figure 7. Comparison of experimental coarsening rate coefficients with model predictions from (a) LSW [126,127], (b) MLSW [132], (c) Ai et al. [137], and (d) Li et al [136].

Equations (11(11) ), (21(21) ), (22(22) ) and (24(24) ) now represent a set of physically based constitutive equations describing the evolution of creep strain rate based on strain-induced mobile dislocation multiplication and thermally induced precipitate coarsening. An amended version of the process flow diagram from Figure 2 which includes these relations is presented in Figure 8. Assumed microstructural parameters are summarised in the inset table. The set of differential equations can be integrated using standard finite difference methods to yield the full shape of the creep curve. This in turn can be compared to experimental results to validate the present model.

Figure 8. Process flow diagram including explicit constitutive equations and microstructural parameters used in the present model.

4. Critical appraisal of model predictions

4.1. Predictions of tertiary creep in CMSX-4

Model predictions are first compared to an experimental creep data set for CMSX-4 [123,142], which includes a wide range of testing conditions over which the model can be assessed. As shown in Figure 9, experimental results are closely reproduced by the present model. Predictions for time to 1% strain – often used as a design parameter for turbine blade applications [143] – are close to measured values. Only two parameters were varied in order to achieve good agreement – the coefficients for dislocation multiplication and attrition, M and φ, in Equation (20(20) ).

Figure 9. Comparison of simulated creep curves with experimental data for CMSX-4 at (a) 850°C, (b) 900°C, (c) 950°C, and (d) 1000°C.

As discussed by Reed et al. [122], the exponential decrease will be more pronounced in the high-temperature rafting regime. This higher value of φ has been related to additional hardening provided by the rafted microstructure, which traps dislocations in the γ matrix channels. In the case of CMSX-4, rafting can be expected to begin playing a role at 950°C and to become dominant at 1000°C [1]. In the present model, φ was taken to be equal to 225 at 850 and 900°C and equal to 1500 at 950 and 1000°C. The former value of 225 corresponds to the attrition parameter obtained by Ram et al. [93] from DDD simulations of tertiary creep in single-crystal superalloy ABD-1 deformed at 900°C and 450 MPa. The latter value of 1500 was chosen based on the attrition parameters derived by Reed et al. [122] from creep tests on CMSX-4 at 1050°C. These figures are assumed to remain constant for the duration of the creep tests for the sake of simplicity. However, the magnitude of operating dislocation attrition mechanisms is expected to fluctuate with changes in alloy microstructure, which will result in varying levels of resistance to dislocation movement [122].

4.2. On the increase in mobile dislocation density

Values of the dislocation multiplication parameter, M, obtained from fitting to experimental data for CMSX-4 are illustrated in Figure 10. Furthermore, a data point obtained from DDD simulation results for ABD-1 deformed at 900°C and 450 MPa has been added [93]. As ABD-1 and CMSX-4 possess similar chemistries and microstructures, comparable rates of dislocation multiplication can be expected. Although the M value from DDD will be strongly influenced by assumptions and restrictions pertaining to this modelling technique, it is nonetheless interesting to note that it is in good agreement with the inferred values for CMSX-4, which gives further confidence in the present approach.

Figure 10. Dislocation multiplication parameters obtained from fitting to experimental data for CMSX-4 plotted against applied stress. A value obtained from DDD simulations of ABD-1 deformed at 900°C and 450 MPa has also been added [93].

Results also confirm the trend reported by Matan et al. [122,142] that the rate of dislocation multiplication increases exponentially with applied stress. Furthermore, the stress exponent in a relation of the form appears to depend on the dominating damage mechanism – in this case, intermediate-temperature dislocation climb at particle interfaces or high-temperature precipitate rafting. In the tertiary regime, a pronounced dependence on applied stress is found at both 850 and 900°C, whereas in the rafting regime at 950 and 1000°C, applied stress plays a lesser role. This is in line with the expectation that the coarsening rate, dictated primarily by temperature, will control deformation rates for rafting.

Quantitative values for and n lack a clear physical meaning, as it has not yet been determined when isolated dislocation sources become active dislocation emitters, nor by which mechanism this occurs. Moreover, previous studies on the influence of stress and temperature on M, and in turn on , in superalloys have only reported the observed trends, but have not discussed these in a physical sense [55,122,142]. For the case of low-temperature plasticity governed by continuous dislocation glide, Kocks et al. [144,145] derived a function of the form: (25) where α is a constant on the order of unity. It is acknowledged that this simplified relation cannot be readily extended to the present case of a complex engineering material deformed at high temperatures. Nonetheless, it suggests a clear dependence of and M on applied stress, as shown in Figure 10. The influence of temperature can be related to a lower shear modulus, G, at higher temperatures.

4.3. Model sensitivity analysis

While model predictions appear reasonable, it is essential to assess their sensitivity to individual parameters. Figure 11 shows simulated creep curves for CMSX-4 at the common design point of 900°C and 400 MPa. In each case, a different parameter is varied while keeping the others constant: (a) volume fraction, , (b) initial mean particle size, , (c) dislocation multiplication parameter, M, and (d) dislocation attrition parameter, φ.

Figure 11. Assessment of the influence of microstructural parameters on the creep performance of CMSX-4 at 900°C and 400 MPa: (a) volume fraction, (b) initial mean particle size, (c) dislocation multiplication parameter, and (d) dislocation attrition parameter.

In Figure 11(a), results for the volume fraction predicted by Thermo-Calc and TCNi8 for CMSX-4 at 900°C, 64.3%, are compared to volume fractions 2% lower and higher. Creep performance is very sensitive to slight changes in fraction. This is in good agreement with the work of Murakumo et al. [146,147] on tie-line single-crystal superalloys. These authors obtained the longest creep lives at 900°C and 392 MPa for 70% alloys and associated this with an increasing Orowan resistance provided by narrower matrix channels, in line with Equation (7(7) ).

The same qualitative trend can be observed in Figure 11(b) for different initial particle size distributions. In this case, the standard deviation of 150 nm was kept constant and only the mean particle size was varied from 450 to 350 and to 550 nm. These distributions were again estimated through histograms containing 5 bins with an equal bin width of 100 nm. As expected, Orowan resistance and creep strength are increased by smaller particles. However, three caveats must be pointed out. First, decreasing initial size becomes increasingly difficult in industrial practice, as it requires higher cooling rates after the solutioning heat treatment, which can result in damaged turbine blades. Second, for most modern superalloys which possess a small, negative lattice mismatch at high temperatures, a decreasing particle size will eventually result in a shape change from cuboidal to spherical particles, which will lower creep strength [58,60]. Third, the importance of initial particle size and of coarsening damage decreases as creep testing conditions become less severe. During isothermal service conditions, all coarsening curves will converge on a certain asymptotic particle size regardless of initial values [27,148]. As such, it can be argued that while it remains important to consider for modelling laboratory creep results, this value is of less importance for industrial purposes.

Figure 11(c,d) reveals an equivalent trend with regard to the evolution of mobile dislocation density, as described by Equation (20(20) ). Creep life is extended as the multiplication parameter, M, increases, or as the attrition parameter, φ, decreases. This observation appears counterintuitive, since a rapid increase in the net rate of mobile dislocations is commonly expected to lead to higher strain rates [13–15]. However, this can be rationalised by considering Equations (21(21) ) and (22(22) ), which describe the net rates of climbing and gliding dislocations, and also the experimental results of Ram et al. [93]. These authors showed that the activation of several individual dislocation sources at low levels of creep strain will result in ‘dislocation avalanches’ within γ matrix channels and, eventually, in the formation of stable, equilibrium dislocation networks. For higher values of M, this drastic surge will be more pronounced and will first lead to a higher initial net rate of gliding dislocations, , in Equation (22(22) ). Second, this also implies a higher initial density of dislocations gliding through the matrix, , which will then become arrested at particle interfaces. Third, following Equation (21(21) ), this will increase the initial net rate of climbing dislocations, , and their density, . Fourth, the higher value of will produce a strong negative feedback effect in Equation (22(22) ) and will outweigh the positive term for the generation of new gliding dislocations, dictated by M. Finally, this effect will markedly reduce the net rate of gliding dislocations, thus inducing an equilibrated state. A lower value of φ will have the same qualitative effect. Overall, this situation can be compared to other physical or biological systems which strive towards a state of equilibrium or homeostasis. In the present work, a higher value of M or lower value of φ will bring about this equilibrium more rapidly by, counterintuitively, causing more damage in very early stages.

Finally, it is critical to examine whether the model can reproduce the increasing creep performance of later-generation superalloys which contain higher amounts of rhenium. The same four compositions presented in Section 3.3.2 are compared at the design point of 900°C and 400 MPa. All microstructural parameters are chosen to be constant, in order to isolate the effects of chemistry on solid solution strengthening and on the rate of particle coarsening. Thermo-Calc and TCNi8 are first used to calculate γ and phase compositions at 900°C. Individual chemistries are then slightly adjusted to yield the same volume fraction of 64.3%, while maintaining the predicted phase compositions.

As shown in Figure 12, creep performance increases as expected in the order SRR99, CMSX-4, TMS-75, and EPM-102. This agrees with experimental results which show significant improvements in creep resistance associated with lower effective diffusivity values and with a retardation of diffusion-controlled deformation mechanisms. However, on the one hand, predicted time to 1% strain is almost identical for SRR99 and CMSX-4, despite the 3 wt.% rhenium additions in the latter. On the other hand, the model underestimates the time to 1% strain for TMS-75 and EPM-102 by about 50%. These discrepancies can be explained by the general inaccuracies associated with CALPHAD calculations and they underline the importance of improving and expanding thermodynamic parameter databases like TCNi8. Nonetheless, alloy performance remains in good qualitative agreement with published data and affirms the usefulness of the present model for alloy design.

Figure 12. Comparison of simulated creep curves at 900°C and 400 MPa for first- to fourth-generation single-crystal superalloys.

5. Discussion: choice of a composition-Dependent creep merit index

One significant goal of the present work is to provide an improved description of creep resistance for use in ABD calculations. In this context, model predictions must be critically assessed for a larger number of compositions and testing conditions. A survey of journal articles, patents, and grey literature yielded around 1500 data points for time to 1% creep strain in different single-crystal superalloys. This set was filtered to include only the tertiary creep regime at intermediate levels of applied stress and temperature, resulting in 265 data points. Upper limits were set at 517 MPa (75 ksi) and 955°C (1750°F), following the results for CMSX-4 shown in Figure 1. Furthermore, alloys predicted by Thermo-Calc and TCNi8 to exceed a total TCP phase volume fraction of 5% were removed, since the effects of these deleterious phases are not considered in the present work. Previous studies showed that a volume fraction of up to 2% TCP phases is tolerable for shorter creep tests with rupture times in the range of 200 to 300 h [149–151], similar to most tests considered here. However, since Thermo-Calc provides estimates for thermodynamic equilibrium conditions, in which a higher volume fraction of TCP phases is expected, the critical limit was increased to 5%.

In Figure 13, experimental results are rendered on the x-axes as temperature-compensated time in the form of the Larson–Miller parameter . Under the assumption that the same damage mechanisms are operating in all the tests considered, the ‘merit indices’ or microstructural parameters rendered on the y-axes should provide relative measures of creep resistance and, in principle, should scale with Larson–Miller values.

Figure 13. Comparison between experimentally measured Larson–Miller parameters for time to 1% creep strain and different ‘merit indices’ for creep strength: (a) effective diffusivity index proposed by Reed et al. [6,56], (b) solid solution strengthening index proposed by Rettig et al. [70], (c) yield stress contribution from solid solution strengthening as calculated by Roth et al. [152], (d) volume fraction predicted by Thermo – Calc, (e) precipitate yield strength as modelled by Crudden et al. [153], and (f) Larson–Miller parameters calculated with the present model.

Previous creep merit indices for alloy design consider a fixed microstructural architecture, with a set volume fraction for all compositions, and neglect coarsening damage. Under this interpretation, creep resistance is primarily dictated by the effective diffusivity controlling climb processes at particle/matrix interfaces. Reed et al. [6,56] proposed the inverse of as a suitable first-order merit index (Figure 13(a)). Similarly, Rettig et al. [70] proposed a weighted solid solution strengthening index which takes into account the relative contributions of molybdenum, rhenium, and tungsten in the γ phase, as derived from experiments on binary Ni–X (X=Mo,Re,W) alloys [59] (Figure 13(b)). When compared to experimental results, neither of these indices appears effective at predicting creep resistance. Rather, data points are clustered together at similar values of the respective index, but are scattered along the x-axes depending on experimental results. A further means of quantifying solid solution strengthening is estimating its contribution towards the total yield stress, as modelled, for example, by Roth et al. [152] (Figure 13(c)). Although this shows somewhat better correlation with experiments, results again appear clustered in a certain range of this index.

The main reason for the poor performance of these three methods is that the underlying assumption of a constant microstructural architecture does not hold true for the range of commercial compositions compared here, or indeed for the range of microstructures one would like to consider during an alloy design process. Both the volume fraction predicted by Thermo-Calc and TCNi8 (Figure 13(d)) and the individual strength of precipitates estimated with the model of Crudden et al. [153] (Figure 13(e)) show large discrepancies within the data set. It can therefore be concluded that a reasonable estimate of creep resistance can only be achieved by considering all contributing factors.

The degree of correlation between experimental results and the factors considered in Figure 13 is illustrated in the scatter plot matrix of Figure 14. To account for the highly non-linear relations, Spearman's ρ rank correlation coefficient is chosen instead of Pearson's r correlation coefficient. Spearman's coefficient assesses the measure of statistical dependence between the rankings of two variables and expresses how well their relationship can be described using a monotonic function. It is less sensitive to strong outliers than Pearson's coefficient, since it limits outlier effects to rank values rather than consider actual numerical values. This is beneficial in the present case, as some strong numerical discrepancies may arise from incorrect CALPHAD predictions. Coefficient values are shown in Table 5. Bold figures indicate correlations significant at the 0.05 level for the experimental results and for the present model.

Figure 14. Scatter plot matrix showing Spearman correlations between experimental results and the six ‘merit indices’ considered previously. The red ellipses indicate 95% confidence levels.

Table 5. Spearman correlation values between experimental results and the six microstructural variables considered previously. Bold coefficient values indicate correlation is significant at the 0.05 level.

	Correlation parameters	Experiment	Creep merit index	Solid solution index	Solid solution strengthening	volume fraction	yield strength	Present model
Experiment	Spearman's ρ	1	−0.38124	−0.04559	0.18384	−0.20658	−0.23577	0.56255
	p-value	–	1.25052E-10	0.45906	0.00261	6.98836E-4	1.03611E-4	1.33641E-23
Creep merit index	Spearman's ρ	−0.38124	1	0.70783	−0.49482	0.64357	−0.06693	0.05797
	p-value	1.25052E-10	–	9.66851E-42	7.79955E-18	1.67699E-32	0.27673	0.34631
Solid solution index	Spearman's ρ	−0.04559	0.70783	1	−0.65094	0.84756	−0.40389	0.43936
	p-value	0.45906	9.66851E-42	–	1.89876E-33	1.327E-74	7.37103E-12	5.57271E-14
Solid solution strengthening	Spearman's ρ	0.18384	−0.49482	−0.65094	1	−0.87384	0.39536	−0.47141
	p-value	0.00261	7.79955E-18	1.89876E-33	–	1.17872E-84	2.19425E-11	4.02182E-16
volume fraction	Spearman's ρ	−0.20658	0.64357	0.84756	−0.87384	1	−0.43714	0.51511
	p-value	6.98836E-4	1.67699E-32	1.327E-74	1.17872E-84	–	7.69633E-14	1.98712E-19
yield strength	Spearman's ρ	−0.23577	−0.06693	−0.40389	0.39536	−0.43714	1	−0.45842
	p-value	1.03611E-4	0.27673	7.37103E-12	2.19425E-11	7.69633E-14	–	3.1566E-15
Present model	Spearman's ρ	0.56255	0.05797	0.43936	−0.47141	0.51511	−0.45842	1
	p-value	1.33641E-23	0.34631	5.57271E-14	4.02182E-16	1.98712E-19	3.1566E-15	–

Statistical analysis confirms that among the six factors considered, the present model offers the highest degree of correlation with experiments, with a Spearman's ρ of 0.56. The present model also shows positive correlations to the solid solution strengthening index proposed by Rettig et al. and to volume fraction. This is unsurprising considering the importance of retarding diffusion-controlled deformation mechanisms and of decreasing the size of matrix channels. Interestingly, no significant correlation can be found for the model and the creep merit index based on effective diffusivity. The negative correlations between the model and the two factors associated with yield strength bring to mind the necessary trade-offs in alloy design studies, in this case balancing mechanical strength at low and high temperatures.

Overall, results from the present model (Figure 13(f)) are shown to be in reasonable agreement with experimental results and largely differ by less than one-half of a Larson–Miller parameter. It is acknowledged that the logarithmic term in the Larson–Miller parameter will blur some differences between model and experiments. Nonetheless, unlike other creep merit indices, the model captures the relative performance of commercial compositions and – from a pragmatic point of view – is a suitable benchmark for Alloys-By-Design calculations. Furthermore, the prediction range of one-half of a Larson–Miller parameter is of the same magnitude as the scatter observed in experimental creep datasets due to (1) variations in chemistry between heats, (2) variations in heat treatment between specimen batches, and (3) variations in crystal orientation between individual specimens [154].

Nonetheless, the influence of several inputs and calculated alloy properties on results is assessed in Figure 15 to rationalise some observed discrepancies. Colour scales illustrate across five evenly spaced major levels whether parameters are beneficial (green) or detrimental (red) to creep life. As the colour maps only offer a qualitative indication of the sign and magnitude of these influences, the following approach is taken to also obtain quantitative values. First, a linear relation is assumed between modelling predictions and experimental results of temperature-compensated time, both expressed in the form of the Larson–Miller parameter . Second, residual values are calculated for each data point based on a linear fit. A histogram of residuals and a normal probability plot of residuals are used to check whether the variance is normally distributed. This is found to be the case, as the resulting histogram shows a symmetric bell-shaped distribution around zero and as the normal probability plot is almost linear. Last, residual values are plotted in Figure 16 against several parameters, in complement to Figure 15.

Figure 15. Qualitative assessment of the influence of creep testing conditions – (a) applied stress and (b) temperature – and calculated material properties – (c) density, (d) volume fraction, (e) coarsening rate parameter, and (f) average metal-d level – on the accuracy of model predictions.

Figure 16. Quantitative residual values, obtained after assuming a linear relation between model predictions and experimental results, plotted against (a) applied stress, (b) density, (c) volume fraction, and (d) coarsening rate parameter.

The present model appears to overestimate the performance of alloys tested at the higher-stress/lower-temperature limit of the tertiary creep regime, as shown in Figure 15(a, b) and Figure 16(a). Test results at applied stresses above 500 MPa have associated residual values of around to units of the Larson–Miller parameter, indicating that model predictions of time to 1% creep strain are too high. At these stress levels, shearing of precipitates becomes feasible [19]. Although this mechanism may not be dominant, it represents an additional source of damage which remains unaccounted for and which explains this overestimation.

Likewise, overpredicted values for volume fraction and/or underpredicted values for the coarsening rate parameter lead to an overestimation of creep life in Figure 15(d, e) and Figure 16(c, d). In both cases, residual values are clearly scaling with these parameters. This highlights the importance of extending and improving CALPHAD thermodynamic databases and coarsening models in order to obtain accurate estimates of microstructural evolution at high temperatures for a wide range of alloy compositions.

As further guidance for alloy design studies, the impact of density and of microstructural stability is considered in Figure 15(c, f). Density is calculated using a corrected rule-of-mixtures approach, which also accounts for bonding characteristics in superalloys [6,56]. For increased performance in service, new alloys must possess higher specific creep strength. This means that the beneficial effects of adding slow-diffusing, but heavy refractory elements like rhenium or tungsten must offset the detrimental effects of higher density and, consequently, of higher stresses arising from centrifugal forces [84]. The average metal-d level, denoted as , provides a basic estimate for the probability of deleterious TCP phase formation. It is calculated with a rule-of-mixtures using the average energy levels of d-orbitals expected for each chemical element from elementary band theory [155,156]. Although it has been established that using a critical cut-off value for as a design criterion is too simplistic [157], this parameter remains a useful first-order estimate of long-term microstructural stability. In an Alloys-By-Design study, compositions combining high specific creep strength and low values would be investigated in more detail. Presuming other design requirements – e.g. oxidation resistance, hot corrosion resistance, castability, and cost – are met, these alloys would become candidates for small-scale trial castings and rapid property assessment campaigns.

6. Summary and conclusions

1. A critical assessment of models for the tertiary creep deformation behaviour of single-crystal superalloys is presented.

2. With the assumptions that the applied stress is too small to enable particle shearing and the temperature is too low to induce directional coarsening and thus the formation of rafted structures, the two main damage mechanisms occurring in the tertiary regime are the coarsening of strengthening precipitates and the accumulation of mobile dislocations in γ matrix channels; creep cavitation does not play a significant role.

3. For accurate modelling, it is shown that it is necessary to account for the influence of mean alloy chemistry on microstructural parameters such as fraction and matrix composition. These serve as inputs for sub-models predicting precipitate coarsening and the accumulation of creep strain.

4. Several methods for calculating coarsening rates in superalloys have been compared. It is shown that using available CALPHAD-based methods, it is possible to accurately estimate both the interfacial energy between the γ and phases and the coarsening rate parameter.

5. Effective drag coefficients for dislocation climb and dislocation glide are demonstrated to be important for prediction of the creep strength imparted by slow-diffusing elements. The ratio of these coefficients is in good agreement with values from literature and can be used to rationalise improvements in creep strength of single-crystal superalloys spanning different generations. These advances are ultimately related to a retardation of diffusional flow at particle interfaces.

6. Good agreement of the present model with the performance of the benchmark alloy CMSX-4 is demonstrated, without the use of excessive fitting.

7. An improved modelling strategy is proposed. One important advantage is the ability to provide reasonable estimates of creep strain life for given applied stress and temperature conditions. A critique using creep data from literature demonstrates that this new approach provides a more accurate merit index for creep resistance, when compared to existing alternatives. It is deemed suitable for alloy optimisation efforts.

Disclosure statement

No potential conflict of interest was reported by the authors.

ORCID

Sabin Sulzer  http://orcid.org/0000-0001-6647-3538

Additional information

Funding

Financial support from Siemens Industrial Turbomachinery AB (Finspång, Sweden) is gratefully acknowledged. RR acknowledges funding from the Engineering and Physical Sciences Research Council (EPSRC) [grant number EP/M005607/01].