A Coupled DEM/SPH Computational Model to Simulate Microstructure Evolution in Ti-6Al-4V Laser Powder Bed Fusion Processes

Abstract

A new multi-stage three-dimensional transient computational model to simulate powder bed fusion (L-PBF) additive manufacturing (AM) processes is presented. The model uses the discrete element method (DEM) for powder flow simulation, an extended smoothed particle hydrodynamics (SPH) for melt pool dynamics and a semi-empirical microstructure evolution strategy to simulate the evolving temperature and microstructure of non-spherical Ti-6Al-4V powder grains undergoing L-PBF. The highly novel use of both DEM and SPH means that varied physics such as collisions between non-spherical powder grains during the coating process and heat transfer, melting, solidification and microstructure evolution during the laser fusion process can be simulated. The new capability is demonstrated by applying a complex representative laser scan pattern to a single-layer Ti-6Al-4V powder bed. It is found that the fast cooling rate primarily leads to a transition between the β and α martensitic phases. A minimal production of the α Widmanstatten phase at the outer edge of the laser is also noted due to an in situ heat treatment effect of the martensitic grains near the laser. This work demonstrates the potential of the coupled DEM/SPH computational model as a realistic tool to investigate the effect of process parameters such as powder morphology, laser scan speed and power characteristics on the Ti-6Al-4V powder bed microstructure.

1. Introduction

Additive manufacturing (AM) based on powder bed fusion (L-PBF) has brought impressive advances in the manufacture of bespoke parts with complex geometries. However, it also poses many technical barriers due to highly transient and varying physical phenomena which occur on a broad range of length and time scales and are difficult to observe and characterize [1]. One of these key challenges is the ability to predict and control the microstructure, and hence the component’s mechanical properties during a L-PBF process. The microstructure is also strongly dependent on thermal history [2] during processing, which, in turn, results from several interdependent and diverse physical processes. This means process optimization and quality control in L-PBF processes is difficult to achieve. Thus, a complete computational model of L-PBF AM that can predict the microstructure evolution and resulting mechanical performance is highly desirable. Before such a model is successfully developed, however, significant hurdles in computational modelling need to be overcome. The difficulties arise from the fact that the AM process is a sum of several sub-processes that occur at different length and time scales and are governed by different physics. These sub-processes include the coating (or spreading) of the powder bed surface, energy transfer from the laser to the powder grains, melting and solidification of the powder, interaction between the melt pool and surrounding solid powder, heat loss to the build plate and microstructure evolution as it cools. In AM, microstructure may be modified during subsequent reheating phases as powder layers above are melted. These processes occur at different length and time scales and therefore require different computational techniques. (It should be noted that in this work the powder particles are referred to as ‘powder grains’ or ‘grains’. The use of the word ‘particle’ is reserved for the basic computational element of the melt pool model.)

Reviews of previous work in modelling various aspects of the powder bed fusion processes (powder coating, melt pool formation and microstructure evolution) are provided in [1,2]. With regard to the powder coating stage, the distribution of powder grains has been simulated using particle packing algorithms [3,4] or the discrete element method (DEM), [5,6,7,8]. Particle packing algorithms place each grain on the previous layer with grains being selectively removed until a desired packing density is reached [3,4]. These methods do not model the inter-granular dynamics during the coating. This contrasts with DEM [5], where the motion and collisions of the grains are governed by inter-granular contact force models. Individual grains have previously been modelled as perfect spheres [9,10,11,12] (which can reasonably represent highly spherical plasma atomized powder) and irregular multi-sphere agglomerates (representing more complex powder shapes such as those found in gas or water atomized powder) [6,7,8,13]. The grains used in our model attempt to replicate the real powder morphology, which has been seen to strongly influence the resultant flow behaviour [14,15,16,17].

The highly complex melt pool formation process has been modelled using traditional mesh-based approaches such as the finite element method (FEM) [18], finite volume methods (FVM) [19] and hybrids of these such as arbitrary Lagrangian–Eulerian (ALE) [20]. Lagrangian (meshless) methods including the lattice Boltzmann method (LBM) [21] and smoothed particle hydrodynamics (SPH) [22] are still relatively novel in melt pool simulations but have significant potential due to their natural ability to represent complex interfacial physics. An extensive review of the current status in melt pool modelling is provided in [23].

Regarding microstructure, various classes of microstructure evolution models exist that can broadly be categorized based on the length scales at which they are applied [24]. Micro-scale models, such as the phase-field method (PFM) [25] can resolve length scales far smaller than the powder grain such as those associated with the growth of secondary dendrites and the dendritic arm spacing. Meso-scale models, such as cellular automata (CA) [26] and kinetic Monte-Carlo methods (MC) [27] are more probabilistic than PFM and typically resolve the growth of the primary columnar crystal structures. Macro-scale models (also known as internal state variable models) [28] operate on the length scale of the thermal field solution. In these models, the microstructure is represented by phase fractions that are computed by solving diffusion type equations that are driven by temperature changes. The representative volume element typically contains many crystals, so the microstructure geometry is not resolved. These models are less precise than the micro-scale and meso-scale models but are computationally more feasible in a fully integrated modelling framework.

While significant progress has been made on modelling the individual physical processes in L-PBF, the development of fully integrated computational L-PBF models is still a key challenge [29]. This is due to the extreme complexity of the L-PBF process which has multiple physical phenomena interacting at spatial and temporal scales that vary significantly. Commercially available computational models of AM focus on the prediction of residual stress and distortion (i.e., the thermal-mechanical response to the laser) and do not consider the melt pool dynamics, the powder coating, geometry of the powder grains or the microstructure formation. Initial work has been done by King et al. [1] in developing a methodology to couple multiple physics-based models on various time and length scales (these include particle packing algorithms for the grain distribution stage and use of ALE for melt pool physics). Chu et al. [30] proposed a concurrent multi-scale method through coupling conservation laws for mass. Molecular dynamics simulations have been coupled to macro-scale FEM [31] and meshfree methods [32] to simulate heat transfer and deformation. In [33] a multi-scale strategy linking a micro-scale phase-field method to a macro-scale FEM was presented and validated for welding applications. These concurrent strategies rely on some type of immediate transfer of information among scales. Choosing what and when to exchange between scales is a key part of these methods.

The main aim of the current work is to present progress towards building a unified model for the L-PBF AM process by linking distinct models of the various sub-processes. At this early stage, while the outputs of some models are used as inputs for other models based on the chronological order of the sub-processes in the L-PBF AM process, the models do not interact with each other during computation, i.e., models are not run concurrently. As part of an effort to build a complete L-PBF AM computational model, a multi-stage coupled computational model using DEM, SPH and a semi-empirical (macro-scale) microstructure evolution model [34] is being developed. The model has been applied to simulate the laser L-PBF process of Ti-6Al-4V powder grains. To model the mechanical coating of the powder bed, a DEM technique is used that can simulate the contact interactions between powder grains of variable size and shape and the interactions with the rake [35], also known as the powder spreader or recoater. In this DEM variant, the powder shape is represented by super-quadrics as this provides the ability to investigate a broad range of granular shapes and aspect ratios with only moderate added computational cost. Grain shape has been recognized as one of the most important parameters influencing the behaviour of granular media; mixing quality and rates, material strength and macroscale flow patterns are all strongly affected by shape [36]. The use of irregularly shaped powders in L-PBF has also been shown to reduce cost [37] and is increasingly being considered for use in commercial AM. Therefore, a broader goal of this work is to understand how powder shape impacts the final products from L-PBF processes. In the current model, satellite particles are, however, not considered.

To model the laser melting of the powder bed grains, SPH is used. SPH was originally developed to simulate astrophysical phenomena [38] but has since been extended to various industrial, geophysical and biological applications [39,40]. While SPH has been applied to solidification processes such as casting [41] its application to L-PBF is still relatively novel. It is a meshless or Lagrangian method in which virtual particles move with the local material velocity and store all fields (such as temperature, velocity and microstructure phases) needed to characterize the material state. This means the heat transfer, phase change and melt flow initiated by the laser passing over the powder bed are naturally modelled with SPH. In the coupled computational model, there are typically hundreds of virtual particles within a single powder grain allowing us to resolve both intra-granular and inter-granular physics for thousands of powder grains in a single simulation. The Lagrangian nature of SPH means that macro-scale microstructure evolution models can be easily applied at each SPH virtual particle. This is exploited in the coupled computational model by applying the microstructure evolution model outlined in [34] at each SPH particle using the particle temperatures as input. This microstructure evolution model is derived from the Johnson–Mehl–Avrami (JMA) theory and predicts transformations between phases as a function of temperature history. This coupling choice means the evolving phase fractions of Ti-6Al-4V can be calculated at each virtual SPH particle. This in turn means the microstructure variation within each powder grain can be represented.

This coupled DEM/SPH computational framework is highly novel and has not been applied to L-PBF applications before. The use of DEM means collisions between non-spherical grains during the coating process can be accurately represented. In addition, the use of SPH means that both intra-granular physics (such as heat transfer, melting, solidification and microstructure evolution) and inter-granular physics (interactions between solid grains and between liquid and solid grains) during the laser fusion can also be represented. A detailed description of the method is first provided and then results presented of the microstructure evolution due to the laser fusion of a non-spherical Ti-6Al-4V powder bed. This work is part of a broader goal of using computational tools to (1) understand/quantify how the microstructure changes with process parameters such as powder morphology, laser scan speed and power and (2) to predict the mechanical performance of the final products.

2. Materials and Methods

In the following sections, the components of the coupled three-dimensional transient computational model are outlined.

Table 1. Definition of symbols and, where applicable, their values, used in this study. Material parameters are obtained from [34,42,43].

Symbol	Description	Units/Values
L	Latent heat release	290 kJ
k	Conductivity	7.2 W/mK
ε	Emissivity	0.32
T	Temperature	K
T0	Initial temperature of build plate and powder bed	300 K
Tsolidus	Solidus temperature	1878 K
Tliquidus	Liquidus temperature	1933 K
Tb	Boiling temperature	3100 K
Tms	Martensite start temperature	848 K
Tref	Build plate reference temperature	300 K
μ	Dynamic viscosity of the melt pool	0.61–10 mPa s
μd	Dynamic friction coefficient of the grains	0.4
en	Coefficient of restitution	0.5
a, b, c	Axes lengths of the super-quadric grains	35–40 μm
m	Shape parameter of the super-quadric grains	2.5–3.0
p	SPH particle separation	5 μm
αms	Ti-6Al-4V martensite α phase fraction	0.0–1.0
αwid	Ti-6Al-4V Widmanstatten α phase fraction	0.0–1.0
αgb	Ti-6Al-4V grain boundary α phase fraction	0.0–1.0
β	Ti-6Al-4V β phase	0.0–1.0
t	Time	s

provides a nomenclature for all symbols introduced in this section.

2.1. Powder Coating (or Spreading) Model

DEM is the ideal tool to simulate the powder coating process, incorporating the dynamics of the powder down to the individual grain level. For this work, the DEM solver developed at CSIRO is used. This solver has been applied to numerous industrial and geophysical processes, see for example [35]. A linear spring and dashpot contact interaction model [44] is used to simulate the collisions between contacting grains. One of the key advantages of this DEM method is its ability to simulate interactions between grains of variable size and non-spherical shape by assuming the grain shape has the following super-quadric form.

$${\left(\frac{x}{a}\right)}^m+{\left(\frac{y}{b}\right)}^m+{\left(\frac{z}{c}\right)}^m=1$$

(1)

Here, a, b, c control the axes lengths of the grain and the exponent m defines the blockiness of the grain (m = 2 corresponds to an ellipsoid shape; as m increases the shape becomes increasingly blocky). For the L-PBF simulations, this means the powder morphology can be accounted for by using a distribution of super-quadric shapes which are matched to the shape distribution of the powder grains of interest. Other material properties such as the friction coefficient and inter-granular cohesion are then calibrated using data from standard laboratory Hall flowmeter and angle of repose tests.

2.2. Melt Pool Model

2.2.1. Extended Smoothed Particle Hydrodynamics (SPH) Method

SPH is used to model the melting of the powder bed due to the laser-induced heating, where the contact region travels over the powder bed. Descriptions of the method and its applications are provided in [38,39,40], so only brief details will be provided here. SPH is a meshless method in which virtual particles (computational nodes) move with the local material velocity. These particles store all fields needed to characterize the material state and its motion. Smoothing kernels (typically spline based) are used to interpolate smoothly varying (differentiable) fields defined everywhere in space from the discrete values of the fields on the surrounding particles. The Navier–Stokes, heat transfer and material response equations are then converted in the SPH discretization process into sets of ordinary differential equations for each field for each particle. These equations are then solved for each particle by explicit second-order integration.

Traditionally, SPH is used to model deforming continuum materials, but in L-PBF applications, the bed consists of solid grains with non-trivial shapes that form a packed bed. Each grain is therefore represented as a discretized cluster of SPH particles. It is only when the grains partially melt that the material is free to flow independently according to the Navier–Stokes equations (as discretized in conventional SPH form). This is sufficient for bed prediction in the melting stage of the process. So, this extended variation of the SPH method has clusters of SPH particles that behave very much like DEM particles, but which are used to solve for heat transfer and can additionally undergo phase changes.

Heat conduction between SPH particles (within grains and between grains at their contact points) follows the method described in [41]. Material surfaces are subject to cooling by radiative heat loss which is calculated for particles at the surfaces of either grains or melt using the Stefan-Boltzmann law with an emissivity ε. The heat transfer at the build plate is modelled using a novel heat conduction sink boundary condition. In this boundary condition, each boundary particle on the build plate is assumed to have a temperature that results from the solution of one-dimensional heat conduction into a semi-infinite body with temperature T_ref (this solution is the complementary error function). The gradient of this temperature function is then used to calculate the conductive heat flux loss at each boundary particle using conductivity k. This enables the size of the computational domain for the melting phase SPH simulation to be minimised whilst having realistic thermal boundary conditions at the edges of the bed.

The heat deposition of the laser is calculated using a ray-tracing method which uses the search grid already defined in the SPH method for efficiently finding neighbour particles, see [38]. The laser is modelled as a set of rays whose density and intensity reflect the laser power and spatial intensity distribution. The intersection between each ray and the first SPH particle encountered is identified. An intersection occurs when the particle position is within a distance p of the closest point on the ray, where p is the SPH particle separation A pre-computed proportion of the associated heat increment is deposited by the ray at that location so the enthalpy and temperature of that SPH particle is consequently increased. It is assumed that all the heat energy from the ray is either reflected from the surface or fully absorbed by that particle since the metal materials are radiatively opaque meaning that the radiative skin thickness is less than the size of one SPH particle.

Recoil pressure from vaporisation of the metal at the laser contact surface is an important fluid dynamic driver of melted metal motion so this is included as a source term in the momentum equation. Its form is dependent on the predicted metal vapour temperature at that location (which is controlled by the temperature of the vaporising metal). Surface tension of the melt and any effects from the flow of the shielding gas are currently not included.

Collisions between SPH discretised powder grains are solved using a similar spring-dashpot contact force model as used in DEM [44]. These contact forces are applied between pairs of SPH particles which are on the surfaces of different powder grains. They are calculated using the dynamic friction coefficient μ_d and the normal coefficient of restitution e_n. It is reasonably assumed that the contacts of non-melted grains are cohesionless. Interactions between a particle in a melted region and a particle on a powder grain surface are calculated assuming that the Navier–Stokes equations are the most suited governing equation for the motion of this liquid phase material.

2.2.2. Melt Pool Material Models

The metal powders considered in this L-PBF study are Ti-6Al-4V alloys and have material parameters whose average values are provided in Table 1. Conductivity, specific heat and viscosity are all temperature dependent with generic forms for the variations assumed. As the material is an alloy, a mushy-zone solidification behaviour is assumed with a latent heat release that varies linearly for temperatures in the range T_solidus < T < T_liquidus. See [45,46,47] for further details of the SPH heat transfer and solidification models and material property variation with temperature. These also include examples of application of this method to melting and freezing of metals during die casting. At temperatures below T_solidus the metal is solid. For higher temperatures, the metal progressively melts. The viscosity variation in the mushy zone is non-linear using the form given in [47]. The timestep for the explicit integration is dependent on the viscosity of the fluid so the dynamic viscosity when the metal is only slightly melted is limited to μ = 10 mPa s. Simulation results are not sensitive to the specific choice as the heating rates in L-PBF process are very high and only a small fraction of material is at these temperatures and only for extremely short periods of time. The viscosity of the fully liquid metal decreases with increasing temperature (according to the data provided in [48]). Since the boiling temperature of Ti-6Al-4V is high (T_b = 3100 K) viscosities are needed for a large range of high temperatures. It is known that the functional dependence of the viscosity at high temperatures has an Arrhenius form. This allows prediction of the decrease from μ = 3.3 mPa s at T = T_liquidus to μ = 0.61 mPa s at the boiling point (T = T_b). The SPH solidification model [46,47] allows individual SPH particles to change state based on their thermodynamic condition. A particle becomes free to move independently of the parent grain as a fluid once the melting condition is met and then freezes in place once its temperature decreases below T_solidus.

2.3. Semi-Empirical Microstructure Evolution Model

The Lagrangian nature of SPH means that it is straightforward to add more physics-based models at each SPH particle. This is a key benefit of the use of SPH for L-PBF modelling. In this work, the SPH melt pool model is combined with a Ti-6Al-4V microstructure evolution model, developed in [34]. This microstructure evolution model represents transformations between phases as either diffusionless (instantaneous) or diffusional (time-dependent) depending on which phases are transforming. The model is a macro-scale (or mean-field) one in that it does not resolve the microstructure crystallography so the effects of heterogenous nucleation, anisotropic crystal growth and impingement are not considered. It estimates the Ti-6Al-4V phase fractions as functions of the current time, temperature and cooling rate. Specific details of the Ti-6Al-4V microstructure evolution model are provided in [34] so a brief description is provided here. The model calculates the volume fractions of the body-centred cubic β phase and three variants of the hexagonal close-packed α phase. If the cooling rate is high, a martensitic α phase is formed (α_ms). For lower cooling rates, a basketweave structure, called Widmanstatten α (α_wid) is created. For even lower cooling rates, another α phase grows at the β grain boundaries (α_gb).

The model assumes transformations between these four phases are either diffusionless or diffusional depending on which two phases are transforming. Figure 1 is a schematic showing the solution process (reproduced from [34]). At each timestep, the current β volume fraction (X_β) is compared against the equilibrium β volume fraction to determine whether to increase or decrease the β fraction. If β decomposes (which occurs during cooling), α_gb is first formed followed by α_wid, as indicated in the left portion of Figure 1. Both these transformations are diffusional. The remaining β transforms to α_ms in a diffusionless process if the temperature, T is lower than the martensite transformation temperature T_ms. If β forms from the existing α phases (which occurs during heating), the dissolution of α_ms to both α_wid and β is first calculated. Then, any α_wid is converted to β and finally, any α_gb is converted. This is indicated in the right portion of Figure 1.

The diffusionless or instantaneous transformations are governed by the Koistinen–Marburger equation [49]. The diffusional (time-dependent) transformations are based on JMA theory. Further details are provided in [34].

3. Results

The potential of the coupled computational model is demonstrated by presenting simulation results of a laser L-PBF process applied to a Ti-6Al-4V powder bed. The temperature history associated with a complex laser pattern applied to a bed that contains a single layer of powder are examined along with resulting microstructural changes.

3.1. Simulation Setup

Figure 2 shows a DEM simulation of the addition of a single layer of powder in an AM machine. A control volume of 1 cm by 3 cm was filled with approximately 700,000 grains under gravity with an initial heap, using periodic boundary conditions across the width of the box. A layer of powder is added by raking a heap across the bed at a velocity of 23 cm/s using a double-toothed recoater. The recoater geometry features two separate rakes spaced 2.5 mm apart, with each rake having 750 μm wide teeth spaced 250 μm apart. The complex dynamics of the powder flow and the ability of the DEM model to capture the detailed interactions with the rake geometry including fine details of the flow of material between the rake teeth can clearly be observed.

The finer details of the melt pool formation and microstructure evolution are considered using a smaller section of powder bed material, created by using a simple ‘deposition from above’ technique to set the initial powder positions and orientations onto the build plate. Using the bed structure from the DEM coating model is straightforward. A plan view and a cross-sectional view of the initial state of the powder bed is shown in Figure 3. The powder bed section is a square tile with side lengths of 1 mm with a reference Ti-6Al-4V metal powder, with material parameters defined in Table 1. The powder grains vary in size and shape with a size range of 35–40 μm diameter. The grains are assumed to be super-quadric shaped with a shape range of 2.5 < m < 3.0 (where m is the super-quadric shape factor defined in Equation (1)). The SPH particles have a uniform separation of p = 5 μm which means each powder grain contains a cluster of approximately 200 SPH particles. The bed has 800 powder grains (276,039 SPH particles) initially.

The initial temperature for the build plate and powder bed is assumed to be ambient (T₀ = 300 K). The initial microstructure composition is assumed to be fully martensitic (α_ms = 1.0). The build plate is assumed to provide a heat conduction sink boundary condition (as discussed in Section 2.2.1) and a no-slip velocity boundary condition.

3.2. Laser Scan Applied to the Powder Bed

The 114 W laser traverses the powder bed with speed in the range 0.6 m/s. The laser traversal takes 14 ms and the total simulation time is 22 ms. The laser has a diameter of 150 μm and a Gaussian intensity profile. The heat energy is deposited to the SPH particles assuming a 35% average absorbance. The laser exposure over the powder bed is a complex sequence of fifteen segments, shown schematically in Figure 4a,b. The first four segments form an outer bounding square and are shown in Figure 4a. The next eleven segments, shown in Figure 4b, create a sawtooth pattern within the region bordered by the outer square. This laser pattern has not been applied in any prior work (either in numerical or physical experiments). The purpose of such an intricate laser pattern is to demonstrate the capability of the coupled DEM/SPH computational model to account for realistic complex laser paths (as used in typical L-PBF operations) and to demonstrate the resulting melting, cooling and remelting in different portions of the powder bed during the laser traversal (which then lead to complexity in the resulting metal microstructure).

3.3. Temperature Evolution of the Powder Bed

The change in temperature of the powder bed as a result of the laser scan is shown in Figure 5a,b. Figure 5 contains eight snap shots of the powder bed temperature during and after the laser exposure. Figure 5a shows the temperature at time t = 1 ms after the first segment of the laser scan (point 1 in Figure 4a). By time t = 4 ms, the laser has traversed the outer square pattern (point 2 in Figure 4a) resulting in the square-shaped temperature distribution shown in Figure 5b. The cooling of the first three segments is now evident at this time. By time t = 6 ms, in Figure 5c, the laser has completed the first six segments (point 3 in Figure 4b). Reheating of the left side of the outer square has occurred as the inner sawtooth pattern has commenced while the other three sides of the outer square have cooled considerably. Figure 5d,e, at times t = 9 ms and t = 12 ms, respectively, show the temperature field at the ninth and twelfth legs of the laser pattern (points 4 and 5 in Figure 4b). The interior of the outer square is progressively heated as the sawtooth pattern continues. The sides of the outer square are now barely evident having cooled back to ambient temperature. By time t = 14 ms, (Figure 5f) the laser has completed its traversal (point 6 in Figure 4b) and the powder bed progressively cools back towards ambient temperature as shown in Figure 5g,h.

3.4. Microstructure Evolution of the Powder Bed

How the evolving temperatures within the powder bed influences the Ti-6Al-4V microstructure is next examined. Specifically the interchange of the β, α_ms and α_wid phases during and after the laser traversal is considered. Note, no formation of the α_gb phase was found to occur in this demonstration. This is likely due to the fast cooling rates in this particular single-layer L-PBF process. (This might be quite different to what would occur during the building of a large part with multiple powder layers and multiple solid layers. Under these circumstances one can expect to see re-heating below the surface and thus the possibility of the formation of α_gb.) Figure 6, Figure 7, Figure 8 and Figure 9 show comparisons of the β, α_ms and α_wid phases at four times during and after the laser exposure. The melted Ti-6Al-4V state (T > T_solidus) is denoted by the grey sections in the β plots.

Figure 6 shows the β, α_ms and α_wid phases at time t = 4 ms after the laser has completed the outer square pattern (point 2 in Figure 4a). Figure 6a–c) shows the β, α_ms and α_wid phases, respectively. Figure 6a,b clearly illustrate the dissolution of the initial α_ms phase to β in grains directly heated by the laser (see the right portion of Figure 1). The top side of the outer square is still hot enough to be in the melted state (T > T_solidus). There is also a small amount of α_wid produced (~0.1% phase fraction) in grains that are at the edge of the laser scan, see Figure 6c. This production of α_wid occurs due to the diffusional heating transformation of α_ms to α_wid (see the right portion of Figure 1). As stated in [34], in this type of transformation, the α_ms actually dissolves to a new morphology with different crystallography and composition to α_wid. For simplicity, however, the microstructure evolution model ignores this new morphology and assumes its crystallography and composition are equivalent to α_wid.

Figure 7 shows the microstructure at time t = 9 ms (point 4 in Figure 4b). At this time, over half of the inner sawtooth pattern has been traversed by the laser. This has resulted in over half of the grains in the interior of the square dissolving from α_ms to β. A significant portion of the square region is still in the liquidus state. We also see evidence of the right side of the initial outer square. This section is yet to be remelted by the laser and has now cooled enough to cause the partial decomposition of β back to α_ms.

Figure 8 shows the microstructure at time t = 14 ms (point 6 in Figure 4b). By this time, the laser has completed its traversal of the bed. The grains in approximately 90% of the interior of the square have completely dissolved from α_ms to β. Cooling of the grains in the outer edges of the square results in the decomposition of β back to α_ms. A distinct outer region of α_wid is also now evident. In Figure 9 at time t = 22 ms, the powder bed has cooled with a significant portion of the grains being converted from β back to α_ms. Little change is noted in the α_wid distribution from this point on as the powder bed cools.

4. Discussion

The key result of this simulation to be explored is the relationship between the processing parameters and the microstructure. Control of the microstructure, and hence the material properties, is the ultimate objective of any coupled model and, of course, of a manufacturer. Thus, the cooling rates extracted from the model and resulting microstructures formed are compared with numerical simulations and experimental observations in the literature.

Powder bed fusion generates cooling rates in the order of 10⁵–10⁶ K/s [50]. The maximum cooling rate, determined both numerically and experimentally, for Ti-6Al-4V depends on the scan pattern and the number of layers printed [50], with a second layer generating lower rates. Cooling rates calculated directly under the laser path of the coupled computational model (compare Figure 5a,b) are of the same order of magnitude, validating our thermal model. At the extreme edge of the laser track the cooling rate is lower, but still significant, in the order of 10⁴ K/s.

Both of these cooling rates are far higher than the critical cooling rate of 410 K/s required to form α_ms in Ti-6Al-4V [51], implying α_ms should be observed throughout the path of the laser after its passage. On transformation from the β phase the microstructures are indeed α_ms in these regions. The occurrence of a very small volume fraction of α_wid (Figure 6c, Figure 7c, Figure 8c and Figure 9c) is the result of the heating of the untransformed α_ms powder to a temperature above T_ms (the martensite transformation temperature) but below the β transus (approximately 1268 K [51]). This phenomenon can be attributed to conduction of heat from nearby grains that are in contact with the laser. This effect is most clearly seen by comparing the heated region of Figure 5a with the same region in Figure 5b. After the passage of the laser, conduction from the powder under the laser path to the adjacent (un-melted) powder results in heating of the adjacent region to 1025 K. This effect is at least consistent with experimental observations of decomposition of martensite in what amounts to an in situ heat treatment-like effect in powder bed fusion [52], except in this case the material in question has not been melted. In practice, the microstructure of un-melted or partially melted powder is not typically examined.

5. Conclusions

A new three-dimensional transient L-PBF computational model that combines DEM, SPH and a JMAK-style mean-field microstructure evolution model to predict microstructure transformation in a L-PBF process has been presented. By coupling these different computational techniques, complex interdependent physical phenomena (such as coating, heat transfer, melting, solidification, deformation of melted metal and collisions of solid non-spherical powder grains) that typically occur in L-PBF processes can be represented. This computational model has been applied to predict the microstructure evolution in a single-layer Ti-6Al-4V non-spherical powder bed during and after the laser scan. While the current aim has primarily been to demonstrate progress towards developing a unified and highly customizable model for the L-PBF AM process by linking distinct models of the various sub-processes, the next stage of the work will involve detailed validation and sensitivity studies of process parameters such as powder morphology, powder bed structure, laser scan speed and power characteristics. When relevant experimental data become available for validation, the parameters used in the current models will be tuned accordingly.

In terms of demonstrating the capabilities, this work has highlighted the fact that the approach taken here is highly amenable to customization. For example, the use of super-quadric shapes for powders shows that the treatments employed can be adapted for individual situations where other probability density functions for different powder shapes are involved. Likewise, the incorporation of well-known models for microstructure evolution points out that other appropriate models may be similarly incorporated for other L-PBF scenarios.