Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences
You have accessResearch articles

A computational study on the three-dimensional printability of precipitate-strengthened nickel-based superalloys

Published:https://doi.org/10.1098/rspa.2018.0295

    Abstract

    This paper presents a computational framework to study the differences in process-induced microvoid and precipitate distributions during selective laser melting (SLM) of two nickel-based superalloys representative of low (IN718) and high (CM247LC) volume fraction precipitate-strengthened alloys. Simulations indicate that CM247LC has a higher propensity to form process-induced microvoids than IN718. Particle sintering is predicted to be strongly influenced by the powder size distribution. For deposition thickness of approximately 40 μm, thermal gradients during cooling are predicted to be larger for CM247LC than IN718 and consequently expect the development of larger residual stresses for a high volume fraction γ′ alloy. A coupled mean field/finite-element approach has been used to predict the precipitate distributions across a simple rectangular build and during a subsequent hot isostatic pressing (HIP) cycle. Unimodal and multi-modal particle distributions are predicted for IN718 and CM247LC at the end of the SLM, respectively. A higher volume fraction of γ′ is predicted for CM247LC at the end of the SLM process. During HIP, simulations indicate a dramatic increase in the γ′ volume fraction in CM247LC, which can result in a reduction in stress relaxation and lead to a ductility drop.

    1. Introduction

    The use of additive manufacture (AM) processes or three-dimensional printing can unlock novel design concepts that would otherwise be considered unfeasible using traditional subtractive manufacturing processes. The aerospace industry, in particular, has a need to use these manufacturing processes in conjunction with materials that are capable of operating in arduous high-temperature environments. The precipitate-strengthened nickel-based superalloys typically used in these applications are traditionally considered as unweldable due to their susceptibility to a variety of mechanisms that result in defects such as strain age cracking, ductility dip cracking, solidification shrinkage and liquation cracking [13]. Unfortunately, metallic powder-based AM processes result in solidification and cooling rates that are far higher than those seen in more traditional welding processes which further exacerbates these defect-generating mechanisms [4].

    Given the inherent flexibility of additive manufacturing processes, there is a potentially huge design space to explore and identify the optimum processing window for avoiding these mechanisms [57]. In addition it is very difficult, if not impossible, to observe the behaviour of the material under representative conditions with appropriate temporal and spatial resolution. As such, it becomes very expensive and time-consuming to do the process optimization experimentally [8] to control microstructure and properties in nickel-based superalloys [912]. However, multi-scale simulation tools offer the potential to further understand the physical mechanisms operating during the AM process. The complex melt flow interactions between heat source and powder-bed have been simulated in order to rationalize the underlying mechanisms responsible for generating defects including porosity [13,14], denudation zone [15], metal vapour micro-jet [16], powder particle dynamics [17] and multi-track defects [18]. Moreover, the modelling of AM-induced grain structure has been studied [1922] to pave the way to explore the design space of potential alloys, manufacturing process parameters, and component geometries. Thus, computational tools are a key enabler for cost-effective, timely entry into service for high-performance additive manufactured parts.

    The aim of this paper is, through computational means, to identify differences in the process-induced microstructure of a low and high volume fraction precipitate strengthened nickel-based superalloys, which can potentially lead to cracking during AM. In particular, what differences in the emerging microstructure in such alloy systems lead to the experimentally observed higher propensity for high volume fraction γ′ alloys to fail. To achieve these aims, a multi-scale materials modelling framework is presented that connects composition and powder size distributions (PSDs) to melt pool physics and the subsequent emerging microvoids (size and orientation distributions and volume fraction) and precipitate size distributions. The numerical study has been carried out on compositions representative of the nickel-based superalloys IN718 and CM247LC. The evolution of the γ′ and γ* precipitates during a subsequent hot isostatic pressing (HIP) cycle has also been studied as a means of gauging the extent of stress relaxation, which may be influential in driving cracking. This paper is structured as follows. Section 2 presents a description of the theoretical framework used for the prediction of process-induced microvoids, which is based on a volume-of-fluid (VoF) approach. The evolution of the precipitate distributions during the deposition process are modelled through a multi-component mean-field approach which is outlined in §3. Numerical implementation of these models and parameters used in the present simulations is described in §4. This is followed by numerical simulations of single-layer depositions in §5 and presents predictions of the microvoid and precipitate dispersion for both alloys under consideration. The evolution of the precipitate distributions are also investigated during subsequent HIP of an SLM build. This is followed by a discussion of the numerical results (§6) and conclusions (§7) from this work.

    2. Modelling solid–liquid–vapour transitions during selective laser melting

    Predicting lack-of-fusion of powder particles or melt-induced trapping of gas/air requires modelling of the heat source/powder interactions and the accompanying solid–liquid–vapour transitions. This section presents a VoF approach that aims to simulate these phase transitions and track interfaces of process-induced voids. Numerical implementation of the approach is covered in §4. Consider a computation domain Ω = Ω1Ω2 that is partitioned into subdomains Ω1 and Ω2 corresponding to the metal solid/liquid and vapour phases, respectively. Let the outer boundary of the computational domain be denoted by ∂Ω(out) and the internal liquid/vapour interface as ∂Ω(int) (figure 1a). At temperatures close to melting, the solid phase will be highly viscoplastic and as such its constitutive behaviour can be approximated as a viscous fluid within the current framework. Thus, the solid and liquid phases within Ω1 are differentiated through the solidus temperature and assignment of appropriate constitutive rules. Suppose that at a point X and time t each phase has densities ρi and velocity fields ui(X, t). Introducing the indicator functions ψi(X, t) for each phase defined as

    ψi(X,t)={1XΩi0XΩi2.1
    and satisfy for all XΩ the condition
    ψ1(X,t)+ψ2(X,t)=1.2.2
    The velocity, u(X, t), and density, ρ(X, t), fields at any point in the computational domain Ω can be expressed as
    u(X,t)=ψ1(X,t)u1(X,t)+ψ2(X,t)u2(X,t)2.3
    and
    ρ(X,t)=ψ1(X,t)ρ1(X,t)+ψ2(X,t)ρ2(X,t).2.4
    From equation (2.1), the gradient of the indicator functions vanish everywhere except at the interface between phases, i.e. the interface ∂Ω(int) is identified by the collection of points for which ∥ψi∥≠0, where the operator is the gradient operator. This is illustrated in figure 1b, which shows a number of lack-of-fusion voids identified by using the ψi(X) functions. The volume occupied by each phase in the computational domain is obtained by integration of the indicator functions, which is also used to determine the volume fraction of process-induced voids. The fields ψi(X, t) and ρ satisfy the following continuity conditions:
    ψ1(X,t)t+.(ψ1(X,t)u(X,t))=m˙vapρ2ψ12.5
    and
    ρ(X,t)t+.(ρ(X,t)u(X,t))=0.2.6
    The sink term on the right-hand side of equation (2.5) represents the phase transition from liquid metal to vapour and m˙vap is given by Geiger et al. [23]
    m˙vap=psatm2πRT,2.7
    where m* is the molar mass, psat is the saturation pressure, R is the universal gas constant and T the temperature. In the present work, it is assumed that psat may be approximated by the recoil pressure precoil using the model proposed by Courtois et al. [24]. Although, each phase is assumed to be incompressible, it is recognized that this condition will break down at the liquid/vapour interface. However, these jump conditions are not taken into accounted in the current study and their impact on the distribution of process-induced lack-of-fusion voids will be addressed in subsequent work.
    Figure 1.

    Figure 1. (a) Definition of computational domain and main physical fields, (b) identification of interfaces using indicator functions ψi.

    The total change in momentum of the system must balance all forces within the continuum. Let the Cauchy stress tensor be σ(X). These will result in surface tractions t(n) on the liquid/vapour interface ∂Ω(int) such that

    (XΩ(int)),t(n)=σn2.8
    where n is the normal vector to ∂Ω(int). Conservation of momentum for the system under consideration reads (refer to figure 1)
    ddtΩρud3X=Ω(outer)t(n)d2X+Ω(int)[[t(n)]]d2X+Ωρbd3X,2.9
    where b is a body forces per unit volume and [[t(n)]] denotes a discontinuity of the surface tractions across the metal liquid/vapour interface resulting from differences in the surface tension across the interface. Let γ be the surface tension on the liquid/vapour interface. It is conventional to express [[t(n)]] in terms of the surface gradient operator acting on the surface tension defined as Gs[γ] = γn(.n) − γ and obtain
    Ω(int)[[t(n)]]d2X=Ω(int)[[σ(n)]].nd2X=ΩGs[γ]ψi(X)d3X.2.10
    Let s be the tangent vector to the liquid/vapour interface. Decomposing Gs[γ] into normal and tangential components
    Gs[γ]=γκn+sγs,2.11
    where κ = .n is the local curvature and ∇s = s.. Combining equations (2.9)–(2.11), the local form of momentum conservation is obtained
    (ρu)t+.(ρuu¯)=ρb+.Sp+(γκn+sγs)ψi(X).2.12
    Denoting the dynamic viscosity of each phase as μi for i = 1, 2 then an effective dynamic viscosity can be defined as μ = ψ1μ1 + ψ2μ2. The deviatoric stress tensor S is related to the velocity field u through the constitutive relation
    S=2μD2.13
    and
    D=12(u+u).2.14
    Contributions to the body force b considered in the present work are: (i) buoyancy driven-flows, b(B); (ii) frictional dissipation force b(f)in the mushy zone. These have the following functional form:
    b(X,t)=b(B)(X,t)+b(f)(X,t)2.15
    b(B)(X,t)=α(T(X,t)Tref)g2.16
    and(T<Ts),b(f)(X,t)=Du(X,t),2.17
    where α is the coefficient of thermal expansion, g is the gravitational acceleration, D is the Darcy drag coefficient, T is the temperature and Tref is a reference temperature.

    The total change of energy of a system is the sum of the internal and kinetic energy balanced by the work done and heat changes occurring within Ω. For the latter, two contributions will be considered. The first of these is associated with the heat input from the laser source rlaser, while the second is the rate of release of specific enthalpy due to vaporization rvap. For simplification of the numerical model, contributions from the work done by the moving interface have not been included in the current model simulations. With these assumptions the global statement of conservation of energy becomes

    ddtΩ(e(X,t)+12ρ(X,t)u(X,t)u(X,t))d3X=Ω((rlaser(X,t)rvap(X,t))ψ1.q)d3X2.18
    from which the local condition is obtained (after inserting Fick's Law for the heat flux q)
    DDt(e(X,t)+12ρ(X,t)u(X,t)u(X,t))=(r(laser)(X,t)r(vap)(X,t))ψ1+k2T(X,t),2.19
    where D/Dt is the material derivative, k is the thermal conductivity and ∇2 is the Laplacian operator. The internal energy is assumed to be dominated by phase transformations and is calculated from the specific heat c¯p(T) as follows:
    e(X,t)=TrefTc¯p(T)dT+h(X),2.20
    where
    h(X)={0T(X)TshmeltingTs<T(X)TvaphvapTvap<T(X),2.21
    where hmelting and hvap are the specific enthalpies of melting and vaporization, Ts is the solidus temperature, Tvap is the vaporization temperature and Tref is a reference temperature. The process of vaporization of liquid metal will absorb heat and an estimate for rvap is given by
    rvap((X,t)=m˙vaphvap.2.22
    The energy input from the laser will be modelled as a volumetric heat source as proposed by Xu et al. [25].

    3. Multi-component precipitation model

    This section describes the framework used in the numerical predictions of precipitation evolution during SLM and subsequent heat treatments. The rapid cooling rates of the SLM process result in a fine distribution of γ′ particles typically approximately 5–10 nm [1]. These dimensions are much smaller than the size of the computational domains considered in the VoF simulations of the heat source/powder interactions outlined in the previous section, which are typically of the order approximately 100–500 μm in the present study. Taking these length scale differences into account, a mean-field description of the particle size distribution is a suitable method for describing precipitate phase evolution during solidification. The mean-field formulation for γ′ precipitation is based on the multi-component approach described by Anderson et al. [26] and further developed in [27] to simulate intermetallic precipitation of IN718 during SLM. The number of particles per volume at a time t with radius r′ lying in the range r ≤ r′ < r + dr is F(r, t) dr, where F(r, t) is the particle radius distribution. The moments of F(r, t) provide information on the precipitate dispersion parameters: the zero moment gives the number of particles per volume Nv, the first moment gives the mean particle size 〈r〉 and the third moment the volume fraction ϕ, i.e.

    Nv(t)=0F(r,t)dr3.1
    r(t)=1Nv(t)0rF(r,t)dr3.2
    andϕ(t)=4π30r3F(r,t)dr.3.3
    The evolution of F(r, t) is governed by the following partial differential equation:
    F(r,t)t+r[F(r,t)r˙(r,t)]=F˙+(r,t)F˙(r,t),3.4
    where F˙+(r,t) and F˙(r,t) are nucleation and dissolution rates. The fast cooling rates following solidification of the melt-pool, fine precipitate particle will nucleate. Because of their small size, minimization of the surface energy will dominate resulting in spherical morphologies. Svodoba et al. [28] proposed a description for the particle growth rate, r˙(r,t), for multi-component systems that takes the functional form [27]
    r˙(r,t)=1r2γS1RTΛ(1rc1r)z(r,t)3.5
    Λ=[i=1m(ckic0i)2c0iD0i]13.6
    rc=2γS2i=1ncki(μkiμ0i)3.7
    andz(r,t)=1+r4πNvr,3.8
    where cki and c0i are the elemental concentrations of species i in the precipitate and matrix, respectively. The D0i are the corresponding matrix diffusivities. γ is the particle surface energy and R is the universal gas constant. The z(r, t) function is a correction term for the growth rate in when dealing with finite volume fractions and is given by Marqusee et al. [29]. The terms S1 and S2 are shape factors which have values of unity when modelling spherical particles, and for cylinders with aspect ratios of h, have the following form;
    S1=S20.881h0.1223.9
    and
    S2=0.2912h2/3+0.5824h1/3,3.10
    where h = H/D, H is the height and D is the diameter of the precipitate.

    The generation rate will be assumed to follow from classical nucleation theory and has the functional form

    F˙+(r,t)=ZβNcexp(ΔGkbT)Pinc,3.11
    where Z is the Zeldovitch factor, β* is the atomic detachment rate, Nc is the number density concentration of nuclei, ΔG* is the Gibbs' free energy for nucleation of a particle, kb is Boltzmann's constant. The incubation probability defined as Pinc=exp(τ/t), where τ refers to the incubation time. To account for complex thermal cycles it is convenient to use the following evolution for the incubation probability [27]
    P˙inc=τteqPinc(1teq+(1ΛdΛdT2rcdrcdT+1σdσdT1)dTdt),3.12
    where teq is the equivalent time to reach the current probability at the temperature of interest and is estimated to be teq =  − ln(Pinc)/τ. The activation energy for formation of spherical particles is
    ΔG=16π3γ3ΔGc2.3.13
    For the atomic detachment parameters the following expression for a multi-component system will be used [28]
    β=4πrc2a4VmΛ,3.14
    where Vm is the molar volume and a the lattice parameter.

    Jou et al. [30] proposed the following Gaussian waveform to approximate the nuclei concentration density,

    Nc=N0Δr2πexp((rrc)22(Δr)2),3.15
    where N0 is the concentration of nuclei sites and Δr is the variance of the nuclei size distribution. An expression for Δr may be obtained from the Zeldovitch parameter. The latter is descriptive of the flatness of ∂G/∂r at r = rc [31]. Δr is approximated as the width of plus or minus a thermal fluctuation kbT from rc,
    Δr=(3Ω2(π)3/21Z)1/3,3.16
    where Ω is the atomic volume. The number concentration of nuclei, N0 is approximated by considering the supersaturation and the mean size of the nuclei (rc),
    N0=ϕn3(ϕeqϕ(t))4πrc3,3.17
    where ϕeq is the equilibrium volume fraction and ϕ is the unit volume containing the particle dispersion. The parameter ϕn describes the fraction of active nucleation sites.

    4. Simulation set-up

    The proposed theoretical framework outlined in §§2 and 3 has been applied to alloy systems IN718 and CM247LC, which have been taken as representative of a low and high γ′ volume fraction nickel-based superalloys, respectively. Their nominal compositions are summarized in table 1. Model parameters used for the heat source model are given in tables 2 and 3. The powder size cumulative distributions for these alloys are shown in figure 2 and are taken from powder suppliers. In this particular case, the IN718 size distribution is more uniform and finer than that for CM247LC.

    Figure 2.

    Figure 2. Powder size cumulative distributions for CM247LC and IN718 used in this study.

    Table 1. Chemical compositions of CM247LC and IN718 (at. %) used in simulations. Nickel is the balancing element.

    alloyAlBCCoCrFeHfMoNbTaTiWZr
    CM247LC5.60.0150.079.38.01.40.53.20.79.50.010
    IN7180.50.0419.018.53.05.10.9

    Table 2. Data used for fluid flow and heat transfer calculations. Note that † is the parameter used in the present work.

    physical propertiesIN718refsCM247LCrefs
    solidus temperature, (K)1533[32]1555[33]
    liquidus temperature, (K)1609[32]1641[33]
    evaporation temperature, (K)31903300
    density of liquid, (kg m−3)7400[32]8250[34]
    atomic mass, (u)57.9459.90
    specific heat of solid, (J kg−1 K−1)625[32]790[35]
    specific heat of liquid, (J kg−1 K−1)725[32]860[35]
    thermal conductivity of solid, (W m−1 K−1)21.3029[34]
    thermal conductivity of liquid, (W m−1 K−1)29.30[32]35[34]
    viscosity, (mPa s)7.5[36,37]8.5[36]
    coefficient of thermal expansion, (K−1)1.63 × 10−5[32]1.82 × 10−5[34,38]
    surface tension coefficient, (N m−1)1.88[39,40]1.82[39]
    temperature coefficient of surface tension,−0.123[39]−1.33[39]
     (mN m−1 K−1)
    enthalpy of solid at melting point, (J kg−1)1.20 × 1061.40 × 106
    enthalpy of liquid at melting point, (J kg−1)1.47 × 1061.70 × 106
    enthalpy change of evaporation, (J kg−1)5.8 × 1055.7 × 105
    atmospheric pressure, (N m−2)101 300101 300

    Table 3. Data used for heat source model in this calculation. Note that † is the parameter used in the present work.

    parametersIN718refsCM247LCrefs
    proportion factor, χ1.01.0
    total beam power, (W)200200
    absorption coefficient (ηlaser)0.26[41]0.26
    beam diameter, (μm)150150
    beam velocity, (mm s−1)30003000

    In this study, simulations of a single powder layer deposition are presented over a computational domain of 250 μm × 1000 μm × 400 μm. The VoF framework presented in §2 was implemented within the C++  open source code OpenFOAM software. Details of the numerical methods used in the OpenFOAM implementation are listed in table 4. To run the precipitation model (outlined in §3) in real time within a finite-element (FE) software (Abaqus), user state variables are used to store details such as the size distribution functions of the particle populations and the matrix and precipitate chemistries. The state variables are initialized to contain the initial precipitate dispersions. A user subroutine has been written to convert the state variables into the required information such as the precipitate size distribution functions, then normalize and reformulate the problem, and solve the continuity method shown in equation (3.4) using a finite difference scheme reported in Anderson et al. [27]. Sub-stepping is used within the user state variable subroutine applying a Courant–Friedrichs–Lewy condition to determine the appropriate time step to evolve the dispersions. For cases where the precipitation details are not needed to determine the material response, the precipitation code is decoupled from the FE analysis, and is applied to simulate precipitation behaviour during the predicted thermal loading as a post-processing operation. The ThermoCalc databases used for thermodynamic and mobility calculations are TTNi8 and MOBNi1, respectively. Mean-field model parameters for both alloys are listed in table 5. For details regarding the shape factor S for IN718 appearing in equation (3.7) refer to Anderson et al. [27].

    Table 4. Numerical schemes used by the VoF implementation within OpenFOAM.

    numerical scheme
    theadend time schemebasic first-order implicit, transient, bounded
    gradient schemesecond-order gradient scheme using face-interpolation
    divergence schemeadvection for kinetic energyfirst order, bounded
    advection of pressurefirst order, bounded
    convectionfirst order, bounded
    momentum fluxsecond order, upwind-biased
    mass fluxlimiter for the vanLeer scheme
    compression termquadratic compression scheme
    turbulence modelsecond order, unbounded
    Laplacian schemesecond order, unbounded, no orthogonal correction
    interpolation schemelinear
    surface normal gradient schemeno orthogonal correction

    Table 5. Mean-field precipitation model parameters.

    IN718CM247LC
    Vm (m3 mol−1)7.30 × 10−61.1 × 10−5
    ϕn(γ′)1 × 10−51 × 10−4
    ϕn(γ*)1 × 10−5n.a.
    ϕn(δ)1 × 10−10n.a.

    Simulations of the manufacture of a rectangular build followed by a HIP cycle have been carried out to study the kinetics and size evolution of precipitates for both alloys. The geometry used in these simulations is shown in figure 3. Thermal and mechanical fields were calculated through a FE implementation of the laser heat source and constitutive behaviour. This was done through a user-defined subroutine (UMAT) implemented in the commercial FE ABAQUS software. The build process was simulated through element activation in the FE model based on the laser position and scanning speed. For the precipitation mean-field model, an in-house FORTRAN code has been developed which is called by the UMAT during the FE analysis. For CM247LC, the HIP cycle simulations assumed a heating rate of 0.05°C s−1 and maintained for 2 h at dwell temperature of 1230°C. A two-step HIPping operation was simulated for IN718, with 0.5 h at 982°C followed by 4 h at 1163°C.

    Figure 3.

    Figure 3. Rectangular geometry used in the simulation of the SLM process.

    5. Numerical results

    (a) Prediction of microvoids

    The results are now presented from simulations focusing on microvoid formation. Figure 4 shows instances of the surface morphology after a laser scan of a 20 μm CM247LC single deposition layer for a 150 W laser power at a speed of 1500 mm s−1, where examples of microvoid formation through lack-of-fusion and melt flow mechanisms are indicated. From these multi-phase fluid dynamic simulations, it is possible to extract information on the microvoid dispersion parameters such as volume fraction, size and shape. Predictions of the microvoid volume fraction as a function of layer thickness for CM247LC and IN718 are shown in figure 5. For each condition investigated, five different instances of the particle deposition were carried out, which resulted in the simulated scatter seen in the microvoid volume fraction. The predicted scatter is smallest for the 20 μm deposition layer and increases with increasing deposition thickness, with the numerical results indicating CM247LC having a wider scatter band in the microvoid volume fraction. In the present study, CM247LC is predicted to have a higher propensity for void formation than IN718 as the layer deposition thickness increases. A typical deposition layers thickness is approximately 20 μm, and it can be seen from figure 5 that CM247LC is predicted to have six times more volume fraction of voids than IN718. This difference is reduced if we consider a deposition layer thickness of approximately 40 μm to three times. The volume fraction of voids is predicted to plateau for both alloys for a deposition thickness greater than 40 μm, at approximately 12 and 4% for CM247LC and IN718, respectively. The plateau in the microvoid volume fraction with increasing deposition thickness is associated with the limited penetration of the metal melt pool depth for the process conditions considered. Figure 6 shows instances of the thermal fields on a cross section of deposited layers for IN718. It can be seen that the simulations predict the melt pool depth does not exceed 40 μm for the conditions investigated. Similar predictions are obtained for CM247LC.

    Figure 4.

    Figure 4. Examples of model prediction of microvoid formation due to (a) lack of fusion of powder particles and (b) melt flow trapping of pores.

    Figure 5.

    Figure 5. Prediction of porosity volume fraction as a function of layer thickness being deposited for IN718 and CM247LC.

    Figure 6.

    Figure 6. IN718 prediction of melt pool depth for as a function of deposition thickness.

    Numerical simulations aimed at assessing the influence of PSD on the void volume fractions have also been carried out. This involved using a finer PSD for CM247LC and a coarser one for IN718. In the present study, this was achieved by swapping the PSDs reported in figure 2. Thus, the CM247LC simulations used the finer IN718 PSD and vice versa. The results are shown as star symbols in figure 5. The effect of using a coarser PSD for the IN718 simulations (red star) is to significantly increase the void volume fraction. In fact, the IN718 predictions for the void volume fraction are comparable to those obtained CM247LC with the same size distribution (black square symbols). A finer PSD result in a reduction in the void volume fraction for the CM247LC parameter set, indicated as the black star in figure 5. Although these results need to be confirmed experimentally, it can be concluded from them that the PSD has a major influence on the formation on the dispersion of process-induced microvoids.

    Simulation of the melt pool and subsequent solidification process allows further information on the morphologies and orientation of voids to be extracted. Such numerical outputs can provide fundamental modelling data for the development of a damage mechanics model for the prediction of cracking conditions at the macro-scale. A useful measure of the void morphology is the sphericity S defined as S = (36πV2)1/3/A, where V and A is the volume and surface area of the void. A void with S approaching unity approximates a sphere. In figure 7, the calculated sphericity distributions for the IN718 and CM247LC for the single layer deposition simulations are shown. Increasing the layer thickness has the effect of spreading the distribution, weighted to higher values of S and increase in the number of microvoids having more equiaxed-to-spherical morphologies. This is more pronounced for IN718 than CM247LC.

    Figure 7.

    Figure 7. Prediction of void shape distribution for IN718 and CM247LC for different thickness of deposition layer.

    The corresponding microvoid size distributions are shown in figure 8. For IN718, a unimodal size distribution is predicted for deposition layer thickness up to 40 μm, but becomes bimodal as the deposition thickness increases beyond this. The predicted void size distribution for CM247LC appears to be more uniform for all deposition thickness's investigated. The development of the bimodal void dispersion in the case of IN718 may be rationalized as follows. After deposition, the powder particles do not pack perfectly resulting in a distribution of voids. The dispersion of voids will be modified following a laser pass due to melting of powder particles and sintering of pre-existing voids. For a given set of process parameters, the melt-pool depth is fixed (figure 6) and when the deposition thickness is large compared to this depth, then only part of the initial packing void distribution is influenced by the heat source. Consequently, a bimodal void size distribution develops: one population (fine void size) is determined by the heat source/powder interactions while the second (coarser-sized voids) are associated with the initial packing. If the deposition layer thickness is smaller or comparable to the heat source length, then it is expected that the resulting void size distribution is dominated by heat/source/powder interactions. An example of the spatial distribution of microvoids obtained from these numerical simulations is shown in figure 9. These correspond to a 40 μm thick deposition layer with void volume fractions of approximately 0.14 for CM247LC simulations and approximately 0.04 for the IN718 predictions.

    Figure 8.

    Figure 8. Prediction of void radius distribution for IN718 and CM247LC for different deposition layer thickness.

    Figure 9.

    Figure 9. Two instances of the predicted microvoid spatial distribution for a 40 μm deposition layer: (a) CM247LC and (b) IN718. (Online version in colour.)

    (b) Precipitate evolution during selective laser melting

    Calibration of the mean-field model parameters for IN718 has been carried out by Anderson et al. [27] and these have been used in predicting precipitation of γ′, γ* and δ in the FE simulations of a rectangular SLM build. The γ* and δ phase precipitates are non-spherically shaped and this has been taken into account in the mean-field model as described in §3. In this paper, we will report predictions for IN718 on γ* only. Calibration of the CM247LC mean field parameters is more challenging due to the lack of experimental data. However, Divya et al. have observed fine γ′ particles with radius of approximately 5 nm. This was the only information used to calibrate the CM247LC mean-field parameters. The mean particle radius and volume fraction are given by the particle size distribution function according to equations (3.2) and (3.3) and these have been calculated as a function of cooling rate. The numerical results are presented in figure 10. The IN718 simulations are consistent with experimental observations [41] with the model predicting both the particle size and volume fraction to decrease with increasing cooling rate.

    Figure 10.

    Figure 10. Simulations of the mean particle size (〈r〉) and volume fraction (ϕ) as a function of cooling rate for γ* in IN718 and γ′ in CM247LC. (Online version in colour.)

    FE simulation of an SLM build has been carried out based on the geometry specified in §4. The predicted thermal fields at three different locations along the axis of build are presented in figure 11. Point p1 is located at the base and point p3 at the top of the build. For the process parameters used in this analysis, the peak temperatures at p1 are approximately 1750°C and approximately 1500°C for IN718 and CM247LC, respectively. The precipitation model predicts a unimodal particle dispersion for the IN718 component, while a multi-modal particle populations is predicted for CM247LC (figure 12). Also, the simulations indicate that the size distribution is location specific, which reflects differences in the thermal history between points of the SLM build. In both cases, larger particles are predicted closer to the substrate, where the component has received longer thermal exposure. The difference is more pronounced in CM247LC compared with IN718, with the largest population near the base having a diameter of 50 nm, while closer to the top the largest population of particles has a mean size of 30 nm. In CM247LC, fine γ′ particles are predicted to have a mean radius of 5 nm. The numerical results for CM247LC are in reasonable agreement with measured values of the γ′ size reported by Divya et al. [1].

    Figure 11.

    Figure 11. Finite-element predictions of the thermal fields at three different location along the build axis. (Online version in colour.)

    Figure 12.

    Figure 12. Simulation of the γ* and γ′ distributions for IN718 and CM247LC, respectively. Predictions correspond to finite-element model in figure 11.

    Changes in mean particle radius and volume fraction of γ* and γ′ in the simulated IN718 and CM247LC builds at location p3 are presented in figure 13. The model predicts nucleation of precipitates upon cooling following a pass of the laser. In some cases, the cyclic nature of the thermal loads associated with multiple passing of the laser, causes dissolution of the precipitated resulting in fluctuations of the volume fraction and particle size.

    Figure 13.

    Figure 13. Precipitation simulations during SLM at point p1 of finite-element model. (a) IN718 γ* evolution and (b) CM247LC γ′ evolution.

    (c) Precipitate evolution during hot isostatic pressing

    The evolution of precipitate distributions for both alloys during a subsequent HIP cycle was investigated. Details of the simulated HIP cycle are outlined in §4. The numerical results of the mean-field model are presented in figure 14. The simulations predict rapid growth and further nucleation of precipitates for both alloys during the initial stages of the HIP loading cycle. However, as the temperature is further increased towards the HIP dwell temperature, dissolution of the precipitates is predicted, with complete dissolution in IN718 while CM247LC still retains a significant amount of γ′. For IN718, the γ* particles reach a maximum volume fraction of approximately 0.13 and a maximum mean radius of approximately 10 nm on thermal loading. The corresponding model prediction for CM247LC indicates that most if not all of the γ′ will precipitate out. This is attributed to the high γ′ solvus temperature for CM247LC. Furthermore, the multi-modal particle populations predicted for the CM247LC build will have slower dissolution kinetics than the unimodal dispersions of IN718, since larger particles have lower growth/dissolution rates than smaller ones (assuming curvature-driven kinetics). A significant growth in the γ′ particle radius is predicted for CM247LC during the dwell stage of the HIP cycle, ranging from approximately 200 to 500 nm.

    Figure 14.

    Figure 14. Precipitation simulations during a HIP cycle at point p1 of the SLM build. (a) IN718 γ* evolution and (b) CM247LC γ′ evolution.

    The flow stress behaviour of nickel-based superalloys is known to be strongly influenced by the presence of second phase particles, since they provide an effective barrier to dislocation movement. The yield stress scales nonlinearly with increasing volume fraction of precipitates. Process-induced residual stresses during SLM provide the necessary driving force for plastic deformation. Consequently, precipitation of particle during the initial stages of the HIP cycle is expected to increase the strength with increasing temperature. For CM247LC, the γ′ never fully dissolves and consequently for this alloy, the relaxation of residual stress will be limited. In comparison, since the precipitates fully dissolve in IN718 by the time the HIP peak temperature is reach, residual stress are allowed to relax further than would otherwise be possible if the particle did not dissolve. Thus, from these simulations, it is anticipated that the combination of a higher microvoid volume fraction and higher strength due to precipitation of γ′ increase the likelihood of cracking in CM247LC compared to IN718.

    6. Discussion

    As already stated, the aim of the present study is, through a physics-based computational approach, provide insights to the differences in the SLM of a low and high volume fraction nickel-based superalloys. In particular, can existing computation methods provide information on the increased propensity of high volume fraction nickel-based superalloys to crack during SLM. It is in this context that we will now discuss the numerical results presented in the previous section.

    The present set of simulations predict a higher volume fraction of microvoids in a CM247LC than IN718 for the process conditions studied. The formation of lack-of-fusion defects is associated with sintering of powder particles neighbouring the melt-pool, which in turn is controlled by the extent of heat transfer away from the liquid metal pool. Changes in energy within the domain are determined by heat transfer mechanisms in the melt-pool and solid state, thermophysical parameters as well as the local packing of the near powder particles. A means of assessing the influence of the local packing of powder particles on the formation of lack-of-fusion microvoids is by changing the PSD. Simulations corresponding to a 20 μm deposition layer were carried out using CM247LC parameters listed in table 2 but with the finer PSD used for IN718, refer to figure 2. Similarly, the corresponding calculation has been carried out for IN718 using the coarser PSD of the CM247LC. The microvoid volume fraction for these simulations are shown in figure 5 and marked by the red and black stars which correspond to the IN718 and CM247LC calculations, respectively. The powder size has a significant impact on the formation of lack-of-fusion and melt-induced defects during SLM. Use a finer PSD for CM247LC has the effect of decreasing the formation of microvoids while increasing the PSD for IN718 is predicted to increase the density of microvoids. Thus, for a given set of process parameters, the minimization of lack-of-fusion voids requires a finer PSD to be used during SLM. However, this will increase the cost of the SLM process, since the atomization of high volume fraction γ′ alloys is expensive.

    Heat transfer within the liquid metal is driven by advective and diffusive transport of energy. A useful measure of the contribution of these mechanisms to heat transfer within the melt-pool is the Peclet number, defined as Pe = L|u|cpρ/k, where L is a characteristic length, |u| the magnitude of the local flow velocity and k is thermal conductivity. Taking L as the melt-pool width, the maximum Pe as a function of deposition thickness has been calculated and plotted in figure 15a. As max(Pe)1 for both alloys, heat transfer in the melt-pool is strongly governed by advective transport. The max(Pe) for IN718 is predicted to be larger than CM247LC, which implies that the IN718 melt-pool has a higher heat content than the corresponding CM247LC case. This is reflected in the simulations of the thermal fields within the melt-pool shown in figure 15b,c for 20 and 40 μm deposition layers. The influence of the PSD on the thermal fields within the melt-pool are shown in figure 15d for the 40 μm deposition simulation. Using a finer PSD for CM247LC and a coarser one for IN718 has the effect of reversing the thermal fields.

    Figure 15.

    Figure 15. (a) Maximum Peclet number as a function of powder layer thickness; (b) temperature profiles within the melt-pool for 20 μm; (c) temperature profiles within the melt-pool for 40 μm; (d) influence of powder size distribution on thermal fields on a 40 μm layer thickness.

    The development of residual stresses during solidification of the melt-pool is governed by thermal gradients. From figure 15b,c, the simulations suggest that increasing the deposition thickness from 20 to 40 μm results in steeper thermal gradients upon cooling from the liquid state for CM247LC than IN718. In practice, deposition layer thickness fall within approximately 30–50 μm and based on these simulation, CM247LC is expected to generate larger residual stresses than IN718. Such stresses can potentially provide the driving force for the development of cracks and/or growth of voids which eventually may lead to failure of the SLM build. The likelihood of cracking is sensitive to the orientation of microvoids relative to the principal residual stress components. Microvoid orientations have been extracted from the present simulations. These are based on fitting an ellipsoid to the process-induced voids and identifying the plane containing the semi-major axis. The normal to this plane is taken to be orientation of the void. The positive direction taken along the normal to the deposition surface The numerical data are presented as a pole figure in figure 16, where the build direction is taken to be [001]. The scanning direction is along the [1¯00] axis on the stereographic projection. From figure 16, it can be seen that CM247LC has a much more randomly orientated distribution of microvoids than IN718 for deposition thickness less than 40 μm. As the thickness of the deposition layer increases the microvoid orientation becomes progressively more anisotropic. This is a reflection of the powder particles arranging themselves more compactly with increasing deposition thickness. The voids will then tend to develop along closed packed planes such as the 〈111〉 orientation.

    Figure 16.

    Figure 16. Pole figure representation of the predicted void orientations for IN718 and CM247LC as a function of deposition layer thickness.

    Mean-field simulations of γ′ precipitation in CM247LC indicate a large thermodynamic driving force for nucleation, resulting in multi-modal particle size distributions. For a low volume fraction γ′ alloy the particle distributions are predicted to be unimodal. a significant amount of intermellaic phase present during the heating up stage of a HIP cycle. This will result in an increased strength of the high γ′ alloy thereby reducing the extent of residual stress relaxation through plastic deformation. For IN718, all precipitate phases dissolve before reaching the peak HIP temperature. Based on the simulations presented this work, it is hypothesized that failure during SLM of high volume fraction γ′ alloys is driven by residual stress acting on microvoids while during HIP of an SLM part, the precipitation of large amounts of γ′ will lead to reductions in ductility. In the present study, it has been demonstrated numerically that a high volume fraction γ′ alloy will be more susceptible to these mechanisms.

    7. Conclusion

    A computational framework for the SLM of nickel-based superalloy components has been presented. The approach couples a multi-phase CFD description of the solid–liquid–vapour transitions accompanying the SLM process to a multi-component mean-field model of the solid-state precipitation reactions.

    Based on commercially available PSDs, simulations of the solid–liquid–vapour transitions for IN718 and CM247LC have been carried out. CM247LC is predicted to have a higher propensity to form microvoids during AM than IN718, with six times more microvoids predicted for a deposition layer of 20 μm. The high volume fraction γ′ alloy is predicted to have larger-sized microvoids.

    The PSD is shown to have a significant impact on the volume fraction of microvoids, with finer distributions resulting in fewer lack-of-fusion and melt-induced voids. The PSD is shown to influence heat transfer within the melt-pool and into the surrounding powder particles. This, in turn, will impact the extent particles can sinter and hence the density of lack-of-fusion/melt-induced microvoids. Within the melt-pool, the heat transfer is shown to be controlled by convection.

    Predicted thermal fields are shown to be influenced by conduction of heat away from the melt-pool. CM247LC has a higher thermal conductivity in the solid state than IN718, resulting in a lower energy content in the melt-pool for a given power input and hence lower temperatures.

    FE simulations of a simple rectangular SLM build were used to predict solid-state precipitation reactions. For CM247LC, multi-modal γ′ distributions were predicted at the end of the SLM process, while a unimodal γ* dispersion was predicted for IN718. Fine particles in the simulated CM247LC build had a mean radius of approximately 8 nm and the large particles were approximately 20 nm. The predicted volume fraction at the end of the build is approximately 0.02. For IN718, a mean γ* particle radius of approximately 1.3 nm is predicted with a volume fraction of approximately 0.07.

    Simulations of the precipitate evolution during a HIP cycle of the simulated rectangular builds have been carried out. For both alloys, further precipitation of phases is predicted at the initial stages of the HIP cycle. At the peak HIP temperature, IN718 all precipitates are predicted to dissolve while in the case of CM247LC a significant volume fraction of particles remain during the HIP dwell.

    Data accessibility

    All data are presented and provided within this paper.

    Authors' contributions

    H.C.B. proposed the modelling framework, managed the work and drafted the manuscript. C.P. performed the CFD simulations using in-house OpenFOAM code whereas Y.S. carried out data processing and analysis using Python scripting. M.J.A. conducted the simulation of precipitation kinetics. R.P.T. supported the finite-element calculations as well as drafting of manuscript. B.S. and J.W.B. coordinated the study. All authors contributed to the manuscript and gave final approval for publication.

    Competing interests

    We declare we have no competing interests.

    Funding

    The authors thank the Aerospace Technology Institute (ATI) for funding this work through the Manufacturing Portfolio programme (project no. 113084).

    Acknowledgements

    The computations described in this paper were performed using the University of Birmingham's BlueBEAR HPC service, which provides a High Performance Computing service to the University's research community. See http://www.birmingham.ac.uk/bear for more details.

    Footnotes

    Published by the Royal Society. All rights reserved.

    References