- Email: [email protected]
A discrete source model of powder bed fusion additive manufacturing thermal history
Accepted Manuscript Title: A Discrete Source Model of Powder Bed Fusion Additive Manufacturing Thermal History Author: Edwin J. Schwalbach Sean P. Donegan Michael G. Chapman Kevin J. Chaput Michael A. Groeber PII: DOI: Reference:
S2214-8604(18)30778-4 https://doi.org/doi:10.1016/j.addma.2018.12.004 ADDMA 603
To appear in: Received date: Revised date: Accepted date:
1 October 2018 7 November 2018 3 December 2018
Please cite this article as: Edwin J. Schwalbach, Sean P. Donegan, Michael G. Chapman, Kevin J. Chaput, Michael A. Groeber, A Discrete Source Model of Powder Bed Fusion Additive Manufacturing Thermal History, (2018), https://doi.org/10.1016/j.addma.2018.12.004 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ip t
A Discrete Source Model of Powder Bed Fusion Additive Manufacturing Thermal History Edwin J. Schwalbacha,∗, Sean P. Donegana , Michael G. Chapmanb , Kevin J. Chaputa , Michael A. Groebera,c,∗ a Air
an
us
cr
Force Research Laboratory, Materials and Manufacturing Directorate, Wright-Patterson Air Force Base, OH 45433 b UES Inc. Beavercreek, OH 45432 c Ohio State University, Department of Integrated Systems Engineering, Columbus, OH 43210
Abstract
Significant attention has been focused on modeling of metallic additive manu-
M
facturing (AM) processes, with the initial aim of predicting local thermal history, and ultimately structure and properties. Existing models range greatly in physical complexity and computational cost, and the implications of various
d
simplifying assumption often go unassessed. In the present work, we first formu-
te
late a fast acting Discrete Source Model (DSM) capable of handling the complex processing often encountered in metal powder bed fusion AM. We then assess
Ac ce p
implications of the source representation, details of the numeric implementation, as well as effects of boundary conditions and thermophysical parameters. We verify the DSM implementation against simple numerical thermal predictions, calibrate it with single track deposit experiments, validate outputs against multitrack deposits, and finally quantify the scaling performance. The DSM is an effective means of quickly generating an estimate of the local thermal history induced by complex scan strategies when combined with arbitrary component geometry. While a number of approximations limit its quantitative accuracy, the inexpensive nature and ability to treat complex processing plans suggests it will be useful for screening and identification of regions experiencing anomalous ∗ Corresponding
author Email addresses: [email protected] (Edwin J. Schwalbach), [email protected] (Michael A. Groeber)
Preprint submitted to Elsevier
November 7, 2018
Page 1 of 48
thermal history. Such a capability is necessary to direct usage of higher fidelity, more expensive models and experimental resources.
ip t
Keywords: selective laser melting, powder bed fusion, thermal modeling,
cr
model verification and validation
1. Introduction
us
Additive Manufacturing (AM) technologies have received considerable research and development attention in the last decade[1, 2, 3, 4, 5]. This collection of materials processing technologies has been touted to provide designers
an
with a number of potential benefits, including manufacturability of unique geometries offering significant improvement in component performance, as well
M
as reductions in either manufacturing cost, lead time, or both. The concept of processing-structure-property relationships is not unique to materials fabricated via AM means, but AM processes can have considerable complexity and offer an
d
extreme degree of location specificity compared to many conventional processes. Therefore, much attention has been focused on modeling of AM processes, with
te
the aim of initially predicting local structure and ultimately properties. With respect to metals AM processes, the local thermal history induced by
Ac ce p
spatially and temporally complex energy input procedures affects microstructure, defects, residual stresses, and ultimately material performance. In this work, we focus on predicting the thermal histories generated by powder bed fusion (PBF) AM processes. Such models have been the focus of intense study over the past several years[6]. Depending on the specific problem, the range of physical phenomena and energy transport mechanisms involved can be extensive: conduction, radiation, phase transformations, and convective transport are all active, with the later driven by capillarity, Marangoni forces, bouyancy, and even recoil pressure[7, 8]. Models incorporating various combinations of such effects have generated considerable insight into the detailed transport phenomena that drive thermal field evolution[9, 10, 11], though high-fidelity solutions to such multiphysics problems are currently still expensive to compute and hence
2
Page 2 of 48
limited in terms of the volume of material they can interrogate. Models that employ reasonable simplifying assumptions in order to efficiently
ip t
incorporate the relevant physics are desirable for screening and trend prediction
purposes, particularly in light of the the widening AM processing parameter
cr
space. Such models inherently rely on some degree of calibration, whether to empirical data or results from higher fidelity models, and must be used with
us
care. It is critical to understand the application intent, range of applicability, and supporting data requirements for such models [12].
The so-called Rosenthal solution is a model originally developed to analyze
an
fusion welding[13], though it has seen much attention from the metal AM community. This is a computationally inexpensive analytic representation of the temperature field induced by a point source of energy moving in a straight line
M
at a constant velocity across a semi-infinite body. Furthermore, it accounts for only conductive transport, and assumes temperature independent thermophysical parameters. As can be inferred from this list of assumptions, this model
d
has several limitations. First, as the energy input is a point source, there is a
te
singularity in the resulting temperature field at the source location. The substrate is taken to be semi-infinite, which may not be unreasonable for very short
Ac ce p
time scales, though it is an increasingly poor assumption as elapsed time becomes long compared the to the thermal diffusion time scale associated with component size. Additionally, the assumption that thermal conduction is the only relevant transport mechanism is questionable in the immediate vicinity of the source and resultant melt pool, where we have already noted that radiation, convective transport, and evaporation can be significant. Finally, and perhaps most importantly, because the energy source moves uniformly in space and time, the model cannot directly describe transient effects induced by source direction changes or modulations in processing conditions such as power or speed, or interactions of multiple sources. In spite of these limitations, this model is useful when operating in the socalled conduction mode. It predicts empirically observed trends in melt pool size and aspect ratio and can be calibrated with an efficiency factor decrementing 3
Page 3 of 48
the input power to bring model output into reasonable quantitative agreement with empirical experience over limited operating parameter ranges. Tang, Pis-
ip t
torius and Beuth derive a melt pool width from the canonical Rosenthal thermal
(1)
us
top surface of the deposit is approximately 1/2 8P w= , eπρcp (Tm − T0 ) v
cr
field[14]. Equation 20 from that work indicates the melt pool width w on the
where P is the absorbed power, ρ is the density, cp is the mass specific heat capacity, Tm is the melting temperature, T0 is the spatially uniform starting or
an
initial temperature of the substrate, and v is the velocity of the beam. Similarly, the length measured from the instantaneous center of the energy source to the trailing edge of the melt pool is found to be
P , 2πκ (Tm − T0 )
M
l=
(2)
where κ is the thermal conductivity. For the case of single track deposits,
te
characterized.
d
all of the parameters appearing in equations 1 and 2 can be reasonably well
In cases where more complex energy input strategies are employed, as is com-
Ac ce p
mon in PBF, the intial temperature T0 can vary significantly with position[15], which in turn has been observed to affect microstructure[16] and residual stress [17]. In principle, the Rosenthal solution could be utilized to estimate thermal field details for cases with more complex energy input strategies, assuming the starting temperature T0 could be easily predicted, and with an additional assumption that the temperature in the local vicinity of the point of interest is sufficiently uniform. However, this is particularly challenging as the energy input path becomes complex, as in the case for typical PBF scan strategies. A number of other researchers have expanded on the model of Rosenthal by replacing the point source with various spatially distributed energy sources [18, 19, 20, 21, 22, 23], effectively regularizing the temperature discontinuity. Additionally, others have focused on expanding the various limitations on path complexity [24, 25, 26, 27, 28, 29, 30, 14, 31]. The main objective in the fol4
Page 4 of 48
lowing work is to not only develop a computationally inexpensive model to predict transient thermal fields produced by arbitrary energy input procedures
ip t
relevant to current and future PBF AM processes, but furthermore to assess model accuracy, establish ranges of applicability, and quantify computational
cr
performance.
The remainder of the paper is organized as follows: We discuss relevant phys-
us
ical assumptions and limitations, formulate a transient thermal history model, and then outline an efficient numeric implementation. We then demonstrate model convergence and validate model output against analytic solutions un-
an
der simple conditions. We calibrate model output against simple experimental results, and then validate model predictions against experimental observations from more complex multipass geometric configurations. Finally, we con-
M
clude with a discussion of model performance and scaling, and then outline future model use in conjunction with reduced order thermal history representations. Throughout this work, we will employ thermophsyical parameters for
d
the aerospace relevant alloy Titanium-6wt%Aluminum-4wt%Vanadium (Ti-6Al-
te
4V). 2. Model Formulation
Ac ce p
2.1. Basic Assumptions
The objective of the present model is to predict an approximate temperature
field during PBF AM processes resulting from (potentially multiple) sources of energy input, such as a laser or electron beam (e-beam), moving along spatially complex paths over a component being fabricated. The primary mechanisms driving thermal field evolution are the energy input itself and resulting thermal conduction.
Energy fluxes due to radiation and latent heat of evaporation can also be significant. Their effects are strongest at very high temperatures and thus their primary influence will be approximately coincident with the molten pool, which is in turn often coincident with the energy input source in both space and time. Laser and e-beam sources can also have absorption efficiencies significantly less 5
Page 5 of 48
than 1.0 [32, 33], and are known to vary appreciably with processing conditions, in particular near the conduction to keyholing mode transition[34]. The present
ip t
model employs an efficiency factor η to decrement the nominal input power P
to account for both the absorption efficiency, as well as losses in the immediate
cr
vicinity of the source due to radiation and evaporation.
Convective flux within the liquid phase is also known to have a significant
us
effect on the shape of melt pool in certain regimes[35, 11]. Rather than directly accounting for these complex physics, the present model modifies the geometry of a volumetric energy source to reproduce the observed melt pool geometries,
an
in particular its depth to width ratio, similar to the approach of Goldak et al.[19]. This shape factor must also be determined via calibration. A process for calibration against experimental data will be discussed in detail in section 4.2,
M
however for both the efficiency and shape parameters, calibration could also be determined by comparison to trusted, higher fidelity models. The latent heat stored or released during liquid-solid and solid-state trans-
d
formations can also affect the thermal field. The net energy contributions from
te
such terms sum to zero over long enough time periods since the amounts of material melted and solidified by the AM PBF process are equivalent. How-
Ac ce p
ever, in the vicinity of the phase transformations these terms will modify the temperature distribution, particularly near the solid-liquid interface. The Stefan number measures the ratio of sensible to latent heat 1 St = Esens /Hsl = Hsl
ZTm cp (T ) dT,
(3)
T0
where Hsl is the mass specific enthalpy of the solid-liquid transformation and cp is the mass specific heat capacity. Heating Ti-6Al-4V from room temperature to melting, we find St = 3.81, using values from [36]. Thus, the sensible heat is roughly four times the latent heat of the transformation, and we expect latent heat to be a secondary, though not completely negligible, effect. We can further compare this to the total amount of energy the beam delivers to the melt pool. Approximating the melt pool volume V as a quarter sphere
6
Page 6 of 48
with radius w/2 and a quarter ellipsoid with axes a = b = w/2, and c = l, where
ip t
w and l are determined from equations 1 and 2, the volume of the melt pool is πw2 w V = +l . (4) 12 2
stored in the liquid is
(5)
us
Esens
ZTm ρ (T ) cp (T ) dT, =V
cr
Assuming the whole volume of this melt pool is at Tm , the total sensible heat
T0
(ignoring the latent heat of solid-state transformations which are typically smaller),
an
and the total latent heat of fusion stored in the melt pool Ef = V ρHsl . Finally, the total energy delivered to the melt pool by the beam is approximately ηP l . v
(6)
M
Eabs =
Substituting typical values for Ti-6Al-4V and assuming an efficiency η = 50% based on reports in references [32, 33, 34], we find Esens /Eabs = 0.275, and
d
Ef /Eabs = 0.071. The total energy delivered to the melt pool is therefore large compared to the energy required to bring it to the melting point, which is in turn
te
large compared to Ef . We will assess the implications of these circumstances further in Section 3.4.
Ac ce p
In light of the above assumptions and simplifications, energy balance in
the presence of N generic volumetric energy sources sˆi which are functions of position ~r and time t, the relevant energy transport equation is N
X ∂ (ρcp T ) = ∇κ · ∇T + sˆi (~r, t) , ∂t i
(7)
where T is temperature, ρ is mass density, cp is mass specific heat capacity, κ is thermal conductivity, and sˆi is power per unit volume introduced by the ith source. With the further assumption that the thermal and transport properties are independent of temperature and uniform in space and time, the above reduces to
N
X sˆi (~r, t) ∂T = α∇2 T + , ∂t ρcp i
(8)
where the thermal diffusivity α is κ/ (ρcp ). 7
Page 7 of 48
2.2. Energy Source Representation
ip t
For conventional PBF AM processes there is a single continuously scanning energy source. To simplify the expression of such a problem, particularly for
cases where the beam moves in arbitrary spatial patterns, we instead represent
cr
this moving source as a series of discrete, stationary energy source terms. The
details of sˆi (~r, t) in equation 8 are then determined by sampling the position of
us
the continuously moving beam, for now at a fixed rate ∆t, and placing stationary sources at these locations ~ri = hxi , yi , zi i, and then activating them at the time τi at which the moving beam passes over them. We refer to this as the Discrete
an
Source Model (DSM), and will demonstrate that on time scales greater than ∆t the thermal profile produced by such a model will converge to that of the original continuously moving source. An example single-track vector is shown
interrogation point of interest.
M
in figure 1(a), with its DSM representation in figure 1(b) near a particular
While formulated in the context of a scanning beam type energy input mode,
d
with appropriate choices of source positions and activation times this method
te
can be used to represent an arbitrary energy input process. It is not limited to a single beam, and in fact could also account for shaped sources such as lines
Ac ce p
or areas. In some cases such as pulsed lasers or the “spot melting” approach employed in [37, 38] to control texture and grain morphology, the discrete approach described here is actually a faithful representation of the energy input process.
The individual source terms employed in the DSM are denoted si and have
a Gaussian distribution in space centered at locations ~ri and supply an impulse of energy at time τi . The effect of a source at location ~ri on the temperature at arbitrary position ~rj and time t is si (Rij , t) =
Ai δ (t − τi ) 3/2
(2πσ 2 )
2 Rij exp − 2 2σ
! .
(9)
Here, Ai δ (t − τi ) is the power of the source with δ the Dirac delta function, σ sets the physical size over which this energy is deposited, and the distance
8
Page 8 of 48
ip t cr 5
10
15
20
25
X[mm]
30
0.025 0.000 −0.025 19.7
19.8
0.0152
(c)
19.9
0.0153
7 5
20.2
20.3
0.0155
0.0156
4
te
Y[mm]
20.1
0.0154
τi[s]
40
d
6
20.0
X[mm]
35
an
(b)
Y[mm]
0
us
1 0 −1
M
Y[mm]
(a)
3 2 1
Ac ce p
0
0
2
4
X[mm]
6
Figure 1: Schematic geometries considered for various problems. (a) 40 mm single vector showing original scan vector (grey arrow) prior to discretization and the interrogation point (black circle). (b) Zoom of the single bead in the region of the interrogation (black circle) showing discrete source locations r~i (colored circles) colored by their timings τi . Open black circles indicate 4σ diameters of two specific sources. (c) 7 × 7 mm square with alternating vectors and several interrogation points.
9
Page 9 of 48
Rij = |~rj − ~ri | is defined for convenience. Note that the units of Ai δ (t) are energy per unit time, and thus the units of si are energy per unit time per
ip t
unit volume. In this model, the source deposits all of it’s energy at t = τi , and
then becomes inactive, though the energy deposited continues to influence the
cr
temperature distribution for all t ≥ τi .
Integrating each si (Rij , t) over the full domain in space and time gives the
us
total amount of energy contributed by the ith source, and equating this total to the effect of the continuous source with efficiency η, τi + ∆t 2
Z Z
Z
ηP (t) dt ≈ ηi Pi ∆t
(10)
an
si (Rij , t) d~rj dt = t
V
τi − ∆t 2
where we have assumed that the continuous source varies in power slowly com-
M
pared to the time scale ∆t and that the effective power of the source ηi Pi is simply equal to the actual power evaluated at τi , ηP (τi ). This approximation is exact when the power is constant over τi − ∆t/2 < t < τi + ∆t/2.
d
For a source located on the object’s top surface defined by zi = 0, we inte-
te
grate over the half space V : −∞ < x < ∞, −∞ < y < ∞, z ≤ 0 and over all time to find the total energy deposited by source i:
Ac ce p
ηi P i
Ai (2πσ 2 )3/2
∆t = R R2 δ (t − τi ) V exp − 2σij2 d~rj dt, −∞
(11)
R∞
Executing the integral and rearranging, we have Ai = 2ηi Pi ∆t.
(12)
The factor of two comes from the fact that half of the total energy introduced by the source is deposited above the plane z = 0. The strength is effectively doubled to ensure that the energy distributed below this plane of symmetry is equivalent to the true continuous source. This also effectively results in a zero-flux boundary condition along the top surface of the part, z = 0.
10
Page 10 of 48
Substituting equation 12 into the original transport equation, we have
+
PN i
= α∇2 T 2 Rij 2ηi Pi ∆t δ (t − τi ) . − 2 3/2 exp 2 2σ ρc (2πσ ) p
(13)
ip t
∂T ∂t
cr
As each new layer in a multilayer deposit is added the coordinate system is adjusted such that the instantaneous top surface of the deposit remains at
us
z = 0. In this reference frame, all sources activated during previous layers will have zi < 0. Such sources can directly influence only material that was added during or prior to the layer during which they became active. Combinations of
an
~ri , ~rj where zi < zj are excluded from the sum in equation 13. This effectively indicates that energy deposited by si cannot diffuse up towards positive z. For cases where the system cools significantly between layers this is not a limiting
M
assumption and the effect of all prior sources is easily captured by the uniform background temperature T0 . However, if the interlayer time is short, or if there is significant non-uniformity of temperature due to heating of newly deposited
d
powder by previously deposited energy, this assumption must be revised. One
te
method to do this would be to accommodate a non-uniform initial temperature condition determined as a function of sources from previous layers, though this
Ac ce p
would require marching the solution in time. 2.3. Thermal Solution
We utilize a Green’s function methodology to obtain solutions to equation
13 in an infinite medium with spatially uniform thermophysical properties and initial temperature T0 [39]. For the case of spherical Gaussian sources we find ! N 2 X Rij −3/2 T (~rj , t) = T0 + Ci Θ(t − τi )λi exp − , (14) 2λi i
where Θ is the Heaviside step function (Θ (x) = 0 if x < 0, and 1 otherwise), and for convenience we have defined Ci =
ηi Pi ∆t √ , ρ cp 2π 3/2
λi = σ 2 + 2α (t − τi ) .
(15) (16)
11
Page 11 of 48
An analogous approach can be taken for a source with an elliptical Gaussian energy distribution. With the ith source dimensions in the x, y, z directions
∆t , 2π 3/2 ρcp
λqi = σqi 2 + 2α (t − τi ) ,
(17)
q = x, y, z,
T (~rj , t) = T0 +
"
N X
(18)
us
the temperature is expressed as
cr
C=√
ip t
respectively being σxi , σyi , and σzi , and defining
Θ(t − τi ) √ ηi Pi C
λxi λyi λzi
i x 2 − 2λijxi
yij 2 2λyi
zij 2 2λzi
an
× exp
−
−
(19)
#
.
where qij = qj − qi for q = x, y, z. In the following work we employ the simplifi-
a plane normal to the beam.
d
2.4. Derivatives
M
cation σx = σy as most laser energy distributions are approximately circular in
~ and cooling rate ∂T /∂t = R are imporThe temperature gradient ∇T = G
te
tant conditions that influence solidification morphology [40, 41]. Solidification 2 ~ , with ~ G front velocity ~v is related to these quantities according to ~v = −RG/
Ac ce p
the sign convention such that the velocity vector is parallel to the thermal gradient in regions that are solidifying (R < 0), and anti-parallel on melting (R > 0). ~ is useful in identifying a position ~r which lies on Additionally, knowledge of G a particular isotherm of interest via Newton’s method. This technique is one method employed to find the position of points on the melting isotherm, and thus measure melt pool dimensions in the following. Analytic expressions for the temporal and spatial derivatives are directly
evaluated from equations 14 and 19, and are listed in the supplemental information for convenience. 2.5. Numeric Implementation Equation 14 indicates that the ith source affects the temperature at all points in space when t ≥ τi . However, the importance of the ith term at an arbitrary 12
Page 12 of 48
location ~rj depends on both the distance Rij and elapsed time t − τi since the source has been activated; increasing either distance or time results in dimin-
ip t
ished influence. Neglecting relatively weak terms in the sum will reduce the computational expense, and, if done with care, will have a minor and measur-
cr
able effect on solution accuracy.
Examining equations 14 and 19, it is clear that when Rij is sufficiently large
us
compared to the denominator of the exponential term, the source’s influence will be weak. Furthermore, σ 2 is typically small compared to the diffusive term as t − τi increases, suggesting that sources which satisfy the simplified condition
an
2 Rij 4α (t − τi ) are relatively unimportant. In practice, a finite cutoff radius
Rcut is chosen, and only terms satisfying the condition Rij ≤ Rcut are retained in the sums in equations 14 and 19. For common scan strategies source terms
M
are approximately uniformly distributed in space, and thus the number of terms retained in the sum will scale approximately as the volume contained within the cutoff distance ∝ Rcut 3 .
d
While the spatial extent of a source term’s influence increases with elapsed
te
time, the magnitude of its impact decreases. Setting the pre-exponential coefficient in equation 14 equal to a minimum temperature change of interest ∆Tmin
Ac ce p
and solving for elapsed time, we have t − τi = Ci / (2α ∆Tmin ). This relationship gives a rational methodology to select a time cutoff tcut , after which source i is no longer important; i is dropped from the sum if t−τi < tcut . We will assess the influence of these approximations on both execution time and solution accuracy in the following sections.
Note that equations 14 and 19 indicate that the computation of T (~rj , t) is
independent of the temperatures at any other points in space or time. Two major implications of this computational decoupling are: 1.) the computation of many space-time temperature outputs is embarrassingly parallel and, as will be seen, scales reasonably well across multiple processors, and 2.) the accuracy of the thermal field approximation is not dependent on the choice space-time output points (~rj , t), i.e. there is no grid stability or output mesh dependence of the solution. 13
Page 13 of 48
ip t cr
us
initialize Tjk = T0 identify active sources for evaluation points (~rj , t) for all evaluation points j [parallel] do for active sources i where Rij < Rcut do for times k where τi ≤ tk < τi + tcut do ∆Tijk = si (Rij , tk ) (equation 9) Tjk = Tjk + ∆Tijk end for end for end for
an
Figure 2: Pseudocode describing the parallel implementation of the present model. Parallelism occurs across the outermost for loop. Index i refers to discrete sources, j to the desired output spatial locations, and k to the desired time output points. Note that the choice of time output points does not have to be equivalent across all spatial output locations.
Parallelization is performed over the spatial evaluation locations, with each ~rj calculation looping only over the source locations ~ri satisfying the time and
M
space cutoffs tcut and Rcut respectively. Each source location evaluates equation 9 for all desired output times of interest. For the specific case of a constant
d
output time interval, the number of time output points is simply tcut divided by the desired output time resolution. Each evaluation of equation 9 is added
te
to a time point for the current spatial evaluation location. Pseudocode detailing this process is shown in figure 2. From the pseudocode,
Ac ce p
it is clear that reducing either the time resolution of evaluations or the source discretization will reduce the computational time of any given spatial evaluation iteration of the parallel loop. An implementation of this algorithm within the DREAM.3D software environment will be described in Section 6[42]. Specifically, a performance study with the source discretization and time resolution are held constant, resulting in each spatial evaluation location in each run of the model having, on average, the same number of evaluations of equation 9. 3. Model Verification In this section, we establish the impact of various physical and numerical approximations and verify the model by reproducing a known analytic result. Section 3.1 assesses the effect of the source discretization, and sections 3.2 and
14
Page 14 of 48
ip t
Description nominal power efficiency beam radius beam velocity melting temp. property temp. mass density heat capacity conductivity starting temp. cutoff distance source discretization cutoff time
cr
Unit W − µm mm s−1 ◦ C ◦ C kg m−3 J kg−1 K−1 W m−1 K−1 ◦ C mm µs s
us
Value 290 1.0 18.75 1300 1650 1203 4252 678.6 22.93 25.0 1.0 20 ∞
an
Symbol P η σ v Tm Tp ρ cp κ T0 Rcut ∆t tcut
M
Table 1: Default parameters utilized for all simulations, unless otherwise noted. The thermophysical parameters are representative of Ti-6Al-4V at a temperature of Tp = 1203 ◦C, as described in [36].
3.6 address numerical implementation effects. Section 3.4 deals with the implications of constant thermophysical parameters and neglect of latent heat.
d
Assumptions regarding boundary conditions and powder bed conductivity are
te
addressed in section 3.5. We defer effects of the neglect of fluid flow until discussing energy source shape in section 4.2.
Ac ce p
Table 1 indicates process, physical, and numeric parameters employed in calculations unless otherwise noted, and figure 1 shows relevant scan geometries that will be considered. Again, thermophysical parameters throughout are representative of Ti-6Al-4V. 3.1. Source Discretization
The relationship expressed in equation 10 guarantees that the same net
energy is deposited independent of the choice of ∆t. However, to accurately
represent the thermal history from a continuously scanning energy source, particularly in the temporal vicinity of an energy input event, an appropriate ∆t must be chosen. Excessively small ∆t will incur increased computational cost and large ∆t will induce spurious thermal excursions. A reasonable guideline to ensure accuracy in the vicinity of an energy input even is that ∆t must be
15
Page 15 of 48
sufficiently small with respect to the time tb = 4σ/v, (where σ = σx = σy in the case of the elliptical source) which is simply the elapsed time required for
ip t
the source to traverse a distance equal to its 4σ beam diameter. It is recognized
that this resolution is higher than needed to represent more spatially distant
cr
input events.
Figure 3 shows the temperature computed at a point on the center of the
us
40 mm long single-bead deposition vector shown in figure 1(a) and (b), as a function of time elapsed since the source was centered on the interrogation point at absolute time t0 . The calculations are performed using the conditions
an
indicated in Table 1, with the exception of the value of ∆t. The solutions exhibit temperature jumps as each discrete source is activated which are strongest in magnitude at times around t0 (see inset). These jumps decrease in magnitude
M
and the solutions become increasingly smooth as ∆t is decreased to a value well below tb . It is also clear that the solutions converge as the elapsed time becomes large compared to ∆t. For the conditions chosen in figure 3, all curves
d
are effectively equivalent approximately 1ms after the beam passes. Unless
te
otherwise stated, for the remainder of this work we employ ∆t = 2 × 10−5 s (or ∆t/tb ≈ 1/3). This discretization is shown in figure 1(b), along with black circles
Ac ce p
representing the spatial extent of two successively activated sources. While there are still visible thermal excursions in figure 3 for the chosen ∆t, the discrepancies
decay at temperatures that are still well above the melting point and thus will not have a significant impact on the predicted melt pool geometry. 3.2. Spatial Cutoff
Next we consider the effect of the spatial cutoff on the quality of the solution.
We employ the same single bead deposit geometry as described in the preceding section, but now vary the spatial cutoff distance Rcut . Figure 4 shows the temperature profile at the center of the bead as a function of time elapsed since the source was centered on the interrogation point at absolute time t0 . The
bounding black curve is the baseline solution with all terms retained in the sum (i.e. Rcut → ∞), and other curves show the result utilizing a range of finite Rcut .
16
Page 16 of 48
105
104 102 −50
0
50 t − t0[μs]
Δt/tb = 3.5 Δt/tb = 0.35 Δt/tb = 0.035 Δt/tb = 0.0035 100
150
103
250
500
750
1000 1250 t − t0[μs]
1500
1750
2000
an
0
ip t
Δt = 2e-04, Δt = 2e-05, Δt = 2e-06, Δt = 2e-07,
103
cr
104
us
T[C]
T[C]
105
d
M
Figure 3: Cooling curves for a point at the center of a single bead deposit using several values of the source discretization time scale ∆t. The beam passes over the point in question at t0 . Inset: expanded time scale surrounding the energy input event. Note: tb = 4σ/v = 58 µs is elapsed required for the beam to traversed it’s 4σ diameter. Nearly complete overlap of conditions 2 × 10−6 s and 2 × 10−7 s indicates convergence with respect to ∆t.
Dashed and dotted curves indicate temperature decay at rates proportional to
te
t−1 and t−3/2 respectively, as a guide to the eye. It is clear that the effect of the spatial cutoff is most siginicant for times
Ac ce p
well after the energy input event. The deviation between the baseline full solution and the various approximations are concentrated primarily at large elapsed times. The accuracy of the thermal history increases as the static cutoff distance is increased. With Rcut = ∞ the solution decays as t−1 for the duration of the calculation (10 s), which is characteristic of diffusive thermal decay in a 2D problem[43]. While the geometry in question is 3D in all cases, the motion of the beam along one direction effectively suppresses thermal diffusion in that dimension resulting in 2D cooling behavior. For cases with finite Rcut , as the elapsed time surpasses a value given by Rcut 2 / (4α) the absence of sources beyond Rcut is “felt” and the temperature decay transitions from cooling proportional to t−1
towards t−3/2 , characteristic of a 3D diffusion case. Ultimately, an appropriate value of Rcut must be chosen to balance compu-
17
Page 17 of 48
105
ip t
104
Rcut = ∞ Rcut = 0.1 [mm] Rcut = 0.3 [mm] Rcut = 1.0 [mm] Rcut = 3.0 [mm] Rcut = 10.0 [mm] ΔT ∝ t −1 ΔT ∝ t −32
101 100 10−1 10−2 10−6
10−5
10−4
10−3
10−2
t − t0 [s]
10−1
100
us
cr
102
101
an
T − T0 [C]
103
M
Figure 4: Cooling curves for a point at the center of a single bead deposit using several values of the spatial cutoff Rcut . The beam passes over the point in question at t0 . Example 3 decay profiles proportional to t−1 (dashed grey) and t− 2 (doted grey) are also indicated for reference.
d
tational expense with accuracy requirements. In the case of figure 4, Rcut =
te
1.0 mm ensures a reasonably accurate thermal history for elapsed times less than 20 ms, at which point the temperature rise induced by the source has dis-
Ac ce p
sipated to approximately 100 ◦C, well below both solid-liquid and solid-state transformations for Ti-6Al-4V. 3.3. Verification: Comparison to Rosenthal Solution The present model can be compared to analytic results for melt pool dimen-
sions in the case of a beam moving at uniform velocity. A series of single bead calculations with power ranging from 150 − 600 W, speed 500 − 2000 mm s−1 ,
and initial temperatures 25 − 1000 ◦C were executed with the DSM. The melt
pool widths were measured by finding points on the melting isotherm when the melt pool was in steady state, far from either the beginning or end of the track, and then compared to results predicted by equation 1. Figure 5(a) shows a comparison of the DSM results with the Rosenthal results and a linear least squares regression between the two. There is excellent agreement with the an-
18
Page 18 of 48
alytic result, with DSM widths generally within 1 µm (.1.0%) of the analytic results across the very wide range of P, v, and T0 considered. A linear regres-
ip t
sion between DSM and Rosenthal widths has slope of 0.996 and an intercept of only 0.35 µm. Similarly, figure 5(b) shows similar excellent agreement for
cr
the melt pool length calculation. Together, these results verify that when run with appropriate numeric parameters the present implementation of the DSM
us
reproduces expected analytic solutions.
For the single track configuration, one other main difference between the DSM and the Rosenthal solution is the finite size of the energy source as char-
an
acterized by the 4σ diameter. We define a P´eclet number for this problem 4σv 4σρcp v = . α κ
Pe =
(20)
M
Employing a 4σ beam diameter on the order of 100 µm and using thermophysical properties for Ti-6Al-4V we find Pe & 10 for v & 1 m s−1 . It is not until the velocity is decreased toward 0.1 m s−1 range that we would expect to see a
d
significant impact of finite spot size on melt pool dimensions.
te
3.4. Thermophysical properties
Another assumption in this model is the temperature independence of the
Ac ce p
parameters κ, ρ, and cp . As indicated in equations 1 and 2, w and l could be
affected by temperature dependence of the product ρcp , and κ respectively. Figure 6 shows the w and l of a melt pool measured at the center of the single bead
deposit described in Figure 1(a) as a function of the temperature Tp at which the fixed parameters are taken for Ti-6Al-4V. For this alloy, the resulting full variation in melt pool widths is on the order of 10% across the full range of property temperature. The length, however, is more sensitive, particularly when using property temperatures below 500 ◦C. This is, of course, alloy dependent, and for the duration of this work we will employ properties at Tp = 1203 ◦C. This value results in effective properties that match the total sensible heat required
19
Page 19 of 48
(a)
cr
700 600 500
us
400 300
1:1 DSM fit: y=0.996x + -0.355 DSM FDM Hsl = 0 FDM Hsl = 286kJ/kg
200 100 0
8000
200
300
400
500
Rosenthal Width [μm]
(b)
7000
600
700
800
d
6000 5000
te
4000 3000 2000
Ac ce p
Numeric Length [μm]
100
M
0
an
Numeric Width [μm]
ip t
800
1:1 DSM DSM fit: y=0.997x + -21.735 FDM Hsl = 0 FDM Hsl = 286kJ/kg
1000
0
0
1000
2000
3000
4000
5000
Rosenthal Length [μm]
6000
7000
8000
Figure 5: Single-bead steady-state melt pool width (a) and length (b) predicted by the DSM and the FDM (with and without Hsl ), vs. the Rosenthal solution[14]. Solutions computed 343 discrete combinations of absorbed power (P = 150 − 600 W), speed (v = 500 − 2000 mm s−1 ), and initial temperature (T 0 = 25 − 1000 ◦C) using thermal properties for Ti-6Al-4V. A linear least-squares regression indicates good agreement between the DSM and Rosenthal solution. The boxed region indicates the region applicable to the experimental results discussed in section 4.1. FDM results computed with Hsl = 0, and 286 kJ kg−1 . Realistic Hsl > 0 has a small impact on steady state width, but a significant impact on the observed melt pool length.
20
Page 20 of 48
3.5
290
3.0
285
2.5
280
2.0
275 270
1.0
Width Length 0
250
500
750
1000 1250 1500 1750 2000
Temperature of fixed parameters [C]
an
265
us
1.5
ip t
295
cr
4.0
Length [mm]
Width [μm]
300
M
Figure 6: Sensitivity of predicted single track melt pool width w and length l as a function of the temperature Tp chosen to evaluate the thermophysical parameters ρ, cp , κ, using parameters for Ti-6Al-4V from Ref. [36]. Vertical dotted lines at T = 1000 ◦C and 1650 ◦C are the approximate locations of the α ↔ β and β ↔ liquid phase transformations respectively.
d
to bring the material from T0 = 25 ◦C to the melting point, i.e.
te
ZTm ρ (Tp ) cp (Tp ) (Tm − T0 ) = ρ (T ) cp (T ) dT.
(21)
T0
Ac ce p
Fully temperature dependent thermophysical properties inherently make the governing equations non-linear, and thus not amenable to the solution approach utilized in the DSM. However, in most alloys this effect is of secondary importance and should not inhibit the use of the DSM for basic initial calculations. We already discussed the effect of Hsl via a simple argument based on the
St number of the alloy. Here, we consider a more complex approach to calculate the effect of the Hsl on aspects of the melt pool shape for the single track case, motivated by results in Section 3.2. It was seen there that for single tracks, diffusion along the direction of beam motion is suppressed, and the temperature initially decays as t−1 , effectlvely a 2D cooling problem. We formulate the diffusion problem in a cylindrical coordinate system (r0 , θ0 , z 0 ), with the cylinder axis z 0 parallel to the direction of beam travel. All derivatives in θ0 and z 0 are
21
Page 21 of 48
Width [µm] 211.1 209.7 211.6 204.8
Width Diff.[%] – -0.7 0.2 -3.0
Length [µm] 1238.9 1210.1 1242.1 1661.7
Length Diff.[%] – -2.3 0.3 34.1
ip t
Configuration Rosenthal DSM FDM, Hsl = 0 FDM, Hsl = 286 kJ kg−1
us
cr
Table 2: Comparison of length and width predicted by the Rosenthal model, the DSM, and the FDM formulated in cylindrical coordinates, in the latter case both with and without latent heat of fusion.
neglected, we set σy = σz , and solve the remaining 1-D diffusion equation with a simple finite difference model (FDM) formulation (see supplemental information
an
for further details). This method is strictly only applicable to sources moving in straight line at uniform velocity, and is included here merely to estimate the magnitude of the effect across a range of vayring P, v, T0 . The apparent heat
M
capacity methodology described in [44] and references therein is employed to account for non-zero Hsl . Table 2 summarizes the predicted steady state w and l of a single track melt pool using the standard conditions with the Rosenthal
d
solution, the DSM approach, and the FDM calculation with both Hsl = 0 as
te
well as a realistic Hsl = 286 kJ kg−1 .
As shown in table 2 and figure 5(a), with Hsl = 0 the FDM solution agrees
Ac ce p
well with the Rosenthal and DSM, the estimated width being generally within 1%, verifying the assumptions and implementation. Using a realistic Hsl = 286 kJ kg−1 in the FDM, the width is slightly decreased and the length and
duration molten grows by ≈ 34% compared to the Rosenthal and DSM lengths, as shown in figure 5(b). Note that in general the effect of Hsl will depend on the specific P, v, and T0 . In particular, neglecting Hsl is a poorer assumption as T0 becomes large because the effective St becomes smaller. For the remainder of the present work, Hsl is ignored. This likely has limited effect on the melt pool width for the alloy in question, but will lead to underestimating the length and melt-duration, on the order of 35%.
22
Page 22 of 48
3.5. Boundary Conditions The solutions in equations 14 and 19 are for a semi-infinite slab with uni-
ip t
form thermophysical parameters throughout the volume. The effective value of α in the powder bed surrounding a component is likely over-estimated in this
cr
configuration, leading to enhanced cooling particularly near the edges of a component. To assess the magnitude of this effect we consider the opposite limiting
us
case by including a zero-flux boundary condition, effectively enforcing zero thermal conductivity in the powder bed outside the part. There are various ways of incorporating boundary conditions into the DSM, but the most expedient for a
an
simple convex cross section is via fictitious “mirror” sources placed outside of the component.
Consider a square cross section 7 × 7 mm as shown in figure 1(c) with con-
M
ventional alternating or “snaking” scan vectors running parallel to x direction and spanning from edge to edge of the square. Figure 7 shows the thermal histories computed at two points on the top surface at positions x, y, z =
d
0.05, 3.50, 0.0 mm (just inside of the boundary), and x, y, z = 2.00, 3.50, 0.0 mm
te
(closer to the center of the feature). The thermal histories are computed with a zero-flux boundary condition applied 75 µm toward the powder-bed side of the
Ac ce p
component-bed boundary. This position is chosen to account for the fact that the melt pool typically extends beyond the position of the center of the beam; some over-melting and finite conductivity in the material immediately adjacent to the component will result. The specific position chosen here is approximately equal to half of a typical melt pool diameter observed for isolated single track cases. Solid lines in the figure are for the no-flux case, and dashed lines are for the original uniform diffusivity case with no additional boundary condition. The difference between the two conditions is shown as a dotted line on the right hand axis. For this analysis, the cutoff distance Rcut is set to a value much larger than the size of the feature such that all locations within the part could be influenced by the boundary. Table 3 summarizes multiple metrics of the difference between histories computed at a number of positions along the center-line of the square. Both the 23
Page 23 of 48
t∆Tmax − t0 [ms] 1.6 2.3 3.2 19.6 53.3 315.8 606.2 722.7
R
∆T dt[◦C s] 248 245 242 236 227 215 202 196
ip t
T∆Tmax [◦C] 1594 1392 1198 662 440 205 151 140
cr
∆Tmax [◦C] 461 341 268 164 105 65 54 51
us
x[mm] 0.00 0.05 0.10 0.25 0.50 1.00 2.00 3.50
an
Table 3: Metrics computed for thermal history differences at several positions x from the boundary for the square geometry described in figure 7, all at y = 3.50 mm). ∆Tmax gives the largest instantaneous temperature difference between the case with and without the no-flux boundary condition. T∆Tmax is the temperature of the no-flux condition solution when the maximum error occurred, t∆Tmax − t0 is elapsed time between the source directly interacting with the point in question and the time when the maximum discrepancy is observed, and R ∆T dt is the integral of the temperature discrepancy over the full history.
M
table and figure 7 indicate that that for a point near the edge, the effect of the boundary condition is generally stronger and is felt at a time closer to the primary energy input event. The maximum difference (∆Tmax = L∞ norm) be-
d
tween the temperature histories is 461 ◦C on the edge of the square (x = 0.0 mm)
te
and decreases toward the center to 51 ◦C. The temperatures T∆Tmax and times t∆Tmax − t0 of the no-flux boundary condition solution at which these maximum
Ac ce p
differences are observed are also tabulated as a function of distance from the edge. It is clear that T∆Tmax decreases as the distance of the position in question from the edge increases, and for distance & 500 µm from the edge, T∆Tmax is
low enough that the effect is likely negligible for this alloy. For the duration of this work we will neglect the effect of the boundary conditions except where stated otherwise, though acknowledge that their influence affects the buildup of heat over long times and their significance is certainly a function of component geometry.
3.6. Time cutoff In Section 3.2 we considered the effect of a varying spatial cutoff Rcut with no time cutoff, tcut = ∞. Here we consider the opposite limiting case with varying tcut and Rcut = ∞. Figure 8 shows the effect of tcut on the thermal
24
Page 24 of 48
ip t cr
us 2500 2000
300
250
T, Temperature [C]
an
T [C]
103
350 T, No flux, x= 0.05[mm] T, No BC, x= 0.05[mm] ΔT, No flux, x = 0.05[mm] T, No flux, x= 2.00[mm] T, No BC, x= 2.00[mm] ΔT, No flux, x = 2.00[mm]
1500
M
1000
500
150 0.150
0.155
0.160
0.165
Time [s]
0.170
0.175
0.180 100
d
102
200
te
50
Ac ce p
0
10−1
100
101
Time [s]
Figure 7: Thermal histories computed at a position near the edge (x, y = 0.05, 3.50 mm, blue) and closer to the center (x, y = 2.00, 3.50 mm, red) of a 7 mm square with a no-flux boundary condition (solid) applied at an offset of 75 µm outside the square, or no boundary condition (dashed) (e.g. spatially constant thermophysical properties in the surrounding bed). Dotted lines indicated the temperature discrepancy, ∆T , between the two solutions on the right hand axis. Inset: zoom of thermal histories in a time window close to the primary energy input event for the positions in question. Note that no distance cutoff is used for this calculation.
25
Page 25 of 48
ΔT, Temperature Error [C]
3000
history for the point at the center of the square feature in figure 1(c), both for the full time history in figure 8(a) and in the region near the primary energy
ip t
input events in figure 8(b). tcut = ∞ is a bounding case with no time cutoff, and the histories increasingly diverge from this case as tcut decreases.
cr
Comparing figure 8 to 4, we can see that finite tcut induces a more rapid error accumulation than Rcut . In the case of finite Rcut , the divergence of the thermal
us
history from the limiting condition at a point ~rj is driven by artificial removal of energy at locations distant from ~rj . This in turn enhances the effective rate of cooling at ~rj , which will, in the limit, lead to cooling that follows t−3/2 behavior.
an
However, in the case of tcut the divergence is a result of artificial energy removal in the immediate vicinity of ~rj , and thus the effect is felt essentially immediately. The rate of error accumulation due to tcut is therefore controlled by the geometry
M
and sequencing of the particular scan strategy employed, not diffusion. Additionally, while the effect of removing an individual source might be small based on the analysis in Section 2.5, if a particular scan strategy results in many
d
points that are spatially and temporally close to each other, their individual
te
weak impacts can still add into large temperature effect. In practice, finite values of both tcut and Rcut are employed, and thus a source could be excluded
Ac ce p
by either threshold. For the reasons noted above, specific choices of tcut and Rcut must be considered jointly, and their their impacts assessed for different scan strategies. In the present case with a single, continuous, and smoothly moving source, the choice tcut = 0.1 s is effective at tracking the history to
temperatures below 500 ◦C, well below the melting point and important solidstate transformations for Ti-6Al-4V. 4. Experimental Measurements Experiments were executed to generate calibration and validation data for the model described above. Parallel single beads were deposited under a range of nine unique P, v combinations, with 2 additional replicates at a centerpoint condition, for a total of 11 beads. An EOS M280 laser PBF machine melted individual vectors after a Ti-6Al-4V powder layer of nominal thickness 30 µm 26
Page 26 of 48
10−1
t [s]
te
1400 1200 1000
Ac ce p
T [C]
100
d
(b)
M
an
102
1600
ip t
us
tcut = ∞ tcut = 1.0 tcut = 0.3 tcut = 0.1 tcut = 0.03
cr
(a)
T [C]
103
tcut = ∞ tcut = 1.0 tcut = 0.3 tcut = 0.1 tcut = 0.03
800 600 400
0.14
0.15
0.16
0.17 t [s]
0.18
0.19
0.20
Figure 8: Effect of varying tcut on the thermal history of the central point in figure 1(c) over the full duration of the scan (a) and a zoom near the most direct energy input event (b). Curves for tcut = ∞, 3.0 s and 1.0 s are effectively equivalent within the window shown in (b).
27
Page 27 of 48
was spread onto a wrought Ti-6Al-4V substrate, initially at room temperature. The individual beads were 40 mm in length and separated by 8 mm. The 4σ
ip t
spot diameter was approximately 75 µm as characterized by beam profiling, and the beam trajectory and speed were verified by an in situ monitoring system[45].
cr
The typical time elapsed between the start of adjacent beads is 0.14 s, and the
total time to deposit all eleven beads was 1.54 s. For comparison, the diffusion
us
time scale l2 / (4α) assuming properties at room temperature for this alloy across a distance l = 8 mm is 2.35 s, and thus the beads are effectively thermally independent on the timescale of the deposition.
an
The resulting single beads were characterized via optical and scanning electron microscopy. Top-down Backscatter Electron (BSE) image montages of single beads were collected in a Tescan Lyra dual-beam FIB. Bead width transients
M
were observed to extend typically no more than ≈ 1mm from the beginning of the bead, and thus the mean and standard deviation of bead width (w, ¯ σw ) were assessed across a steady-state region of interest between 2 and 8 mm from the
d
beginning of the beads. The width was recorded at approximately 30 discrete
te
locations each separated by 200 µm.
Additionally, approximately 15 cross-sections were collected within the same
Ac ce p
steady-state region for each bead via serial sectioning [46]. Cross-sections were etched with Kroll’s reagent (92% H2 O, 6% nitric acid, 2% hydroflouric acid) and optically imaged to reveal the melt pool and surrounding heat affected zone. The mean and standard deviation of the width (w, ¯ σw ), depth of the ¯ σd ), and height above minimum melted location below the substrate surface (d, ¯ σh ) were characterized from these images by three independent the substrate (h, observers, and the results were averaged together. Several examples of crosssection images are shown in the panels of figure 10. The typical variation in the mean width as determined via the top-down compared to cross sections was approximately 1%, and for no bead was the variation > 3.5%. Figure 9 shows a comparison of the top down and cross section measurements, along with bars indicating the measurement standard deviations.
28
Page 28 of 48
180
y
= 0.966x + 3.231
ip t
160
cr
150 140 130 120
130
140
150
160
Cross Sectional Width [μm]
170
180
an
120
us
Top Down Width [μm]
170
M
Figure 9: Bead width as measured via cross sections and top down views. Bars indicate 1 standard deviation on either side of the mean value.
4.1. Single Track Width and Depth
d
Experimentally measured width and depth are shown in Table 4. As ex-
te
pected, both width and depth increase as either power is increased or speed is decreased. Additionally, the depth to width ratio ranges from 0.678 to 0.921 (a
Ac ce p
circular cross section is characterized by d/w = 0.5). In general, the depth to width ratio increases with increasing power, or decreasing speed. The value of the source depth to radius utilized in the model, σz /σx , must be calibrated to
match the observed depth to width ratio. The choice σz /σx = 1 results in melt
pools with circular cross sections and thus σz /σx must be increased to match
the observed melt pool cross section behavior. Furthermore, experimental data indicates considerable variation in w at a
fixed P/v. In particular, 3 replicate tracks at the center-point P, v conditions were produced, and observed mean widths are shown in Table 4, as condition IDs 1, 2, and 9, annotated with asterisks. Tracks 1 and 2 demonstrated similar widths, whereas track 9 exhibited a width ≈ 10% lower. For the purposes of calibrating the model and computing error between empirical observations and model outputs, we utilize the values from track ID 2 for this particular P, v 29
Page 29 of 48
ip t cr us an M d te Ac ce p
Figure 10: Examples of single bead cross sections utilized to determine width and depth. (a) condition 11 (218 W, 975 mm s−1 ), (b) condition 2 (290 W, 1300 mm s−1 ), and (c) condition 10 (363 W, 1625 mm s−1 ). Approximately 30 top down measurements and 15 cross section measurements were collected for each P, v condition.
30
Page 30 of 48
Depth[µm] d¯ σd 100.4 9.2 112.1 4.3 95.3 3.1 141.5 1.8 90.0 4.2 133.9 5.8 86.2 3.7 150.2 3.9 116.9 2.8 121.2 3.7 109.3 3.1
Height[µm] ¯ h σh 25.9 5.7 21.7 3.4 29.4 24.0 20.7 3.2 19.0 6.7 22.1 1.9 26.4 9.8 16.0 4.0 13.6 4.6 22.1 11.1 11.5 2.0
ip t
Width: Top[µm] w ¯ σw 144.1 6.5 143.3 7.1 125.4 8.4 155.4 6.3 125.4 5.3 175.2 5.4 124.4 6.4 172.0 6.4 129.6 4.8 130.6 5.5 152.0 3.8
cr
Width: Sec.[µm] w ¯ σw 147.9 7.0 144.9 5.2 127.2 5.6 153.9 4.2 126.0 5.4 177.3 6.7 122.2 4.8 174.2 8.0 134.0 4.8 131.6 3.5 156.2 4.7
v[mm s−1 ] 1300 1300 1860 1300 1300 908 1733 1040 1300 1625 975
us
P [W] 290 290 341 341 239 239 290 290 290 363 218
an
ID 1* 2* 3 4 5 6 7 8 9* 10 11
M
Table 4: Mean and standard deviation across all measurements of the width (w, ¯ σw ), depth ¯ σd , and height h, ¯ σh of single tracks 1 through 11, in micrometers. Width is characterd, ized by via both cross section (Sec.), and top-down images (Top), whereas height and depth are determined from cross section images only. * denotes replicates at the center-point of the processing conditions.
combination, as it is the most regular of the three deposits. Further work is
d
needed to understand the source of bead-to-bead variation in cases with nominally identical P and v values, but we note that these differences are still small
te
compared to the overall range of widths observed across the full range of P and v in the experiments.
Ac ce p
4.2. Model Calibration
As described in Sec. 2.1, two adjustable parameters are available to bring
the predicted melt pool geometry into agreement with observations. Specifically, the overall efficiency factor η and the volumetric source’s depth to radius ratio
σz /σx must be adjusted to reproduce single track width and depth observations as described in table 4.
Calculations of the melt pool width and depth were performed with the
DSM using the same geometry and processing conditions as the experimental single beads described above. The melt-pool dimensions are determined from extents of the isotherm at Tm as in section 3.3. A range of η and σz /σx values were employed, and width predictions resulting from optimal choices of these calibration parameters are presented in figure 11 along with experimentally 31
Page 31 of 48
Experimental
180
cr
Width
[μm]
200
σz/σx only η only σz/σx and η
ip t
220
160
120
0.18
0.20
0.22 0.24 P/v [J/mm]
0.26
0.28
0.30
an
0.16
us
140
M
Figure 11: Single bead measurements and calculations with varying degrees of calibration. Linear fits over a narrow range of P/v are included as a guide to the eye.
observed widths as a function of linear energy density P/v. Note that P is the
d
nominal power setting before any efficiency factor is applied.
te
Figure 11 shows the model performance under three different calibration cases. In the first case σz /σx is adjusted at constant η = 1 until the depth/width
Ac ce p
ratio best matches observations at the center-point in the P, v space. In the second case, η is tuned to match the observed width with fixed σz /σx = 1.0. In the final case, σz /σx is first adjusted until the depth to width ratio matches observations at the center-point in the P, v space and then η is adjusted until the width matches across the P/v range. The best agreement between modeled and observed width occurs when both
η and σz /σx are tuned, as indicated in figure 12. It should be noted that the optimal choice of σz /σx varies as a function of P and v, but here we only optimize for the centerpoint of the P, v values. Note that equation 1 indicates that w scales as
p P/v. However, the range of
P/v in figure 11 is relatively narrow and does not exhibit extensive curvature and thus linear fits are used to characterize the relationship. Care must be exercised if extrapolating beyond the range of observed data, though calculations with 32
Page 32 of 48
190
ip t
180
160
140
cr
150
us
/
σz σx only η only σz σx and η
130
/ 1:1 1:1±5,10%
120 110 110
120
130
140
150
160
Experimental Width [μm]
170
180
an
Model Width [μm]
170
M
Figure 12: Experimental vs. modeled widths for single tracks across all power, speed conditions for several values the efficiency factor η, and source depth to radius ratio σz /σx .
d
the DSM over a much wider range of P and v values have been executed and p confirm the expected P/v scaling.
te
As expected, the model over-predicts the bead width significantly with η = 1.0. Decreasing the efficiency while keeping σz = σx brings the predictions in
Ac ce p
reasonably good agreement with the observed widths, though the slope of the width with respect to P/v is underestimated. Finally, tuning both η and the σz /σx ratio to match the width at the center-point and brings the slope of the w vs. P/v into better agreement with observations. Again, we point out that previous experimental work has indicated that
the absorption efficiency can be a strong function of the processing conditions, particularly as power is increased near the onset of keyhole instability [34]. Thus, we do not expect a single combination of η, σz /σx to work over a wide range of processing conditions. Instead, the model must be calibrated over the full intended use range of P and v. For the remainder of this work, we will
employ P and v conditions representative of the centerpoint in the range shown in figures 11 and 12, and thus will utilize the η and σz /σx tuned to the center point of this range. 33
Page 33 of 48
5. Model Validation
ip t
After calibrating the DSM against single track experimental results, it is tested against empirical observations of melt pool geometry on the top surface
of an object produced by a multi-vector scan pattern. The test geometry is a
cr
rectangular prism with dimensions of 40 × 10 mm in the build plane, and 20 mm
along the build direction, shown in figure 13. The scan pattern used to process
us
the topmost layer of the item is shown in figure 13(a), colored by the position of the laser as a function of time. The top surface is divided into sub-domains known as ‘stripes’ of width 7 mm. The stripe boundaries are inclined to the long
an
axis of the component by 68.26◦ , and are clearly visible in the BSE montage of the object as a series of dark lines crossing the item in figure 13(b). The
M
individual scan vectors within the boxed region of interest are orthogonal to the stripe boundaries, and are shown in more detail in figure 14(a). When the beam reaches the end of each vector it is shut off before turning and accelerating to
d
begin scanning the next vector, at which point it is turned back on. The mean empirically observed turning time is employed for all turning events in the path,
te
and the beam is assumed to move at the nominal processing speed as it traverses from the final vector in a stripe to the first vector in the subsequent stripe.
Ac ce p
A calculation with this scan pattern is executed with the parameters listed
in table 1 with the exception of the following modifications; the spatial cutoff is set to Rcut = 4 mm to ensure the boundary condition can affect internal regions consistent with sections 3.2 and 3.5, tcut = 0.1 s as described in section 3.6, and T0 = 100 ◦C which is representative of a multilayer deposit. A no-flux boundary condition is employed at a distance 150 µm outside all free edges of the item. The specific region of interest shown in figure 14 is selected because the scan-
ning vectors interact with both a free surface (the lower edge), as well as the internal stripe boundary, annotated as a solid black line labeled “0”. The vectors are sequentially scanned beginning in the upper right of the region and moving toward the lower left. While the processing velocity is nominally uniform, because the vectors become successively shorter as the corner is approached the
34
Page 34 of 48
ip t cr us an
M
Figure 13: (a) Global scan path covering a 40 × 10 mm rectangular prism with scan vectors inclined 21.74◦ from the horizontal. Color indicates the time elapsed since the beginning of the scan on the layer. (b) BSE image montage covering the full object. Boxes in (a) and (b) indicate the region of interest shown in Figure 14, where individual scan vectors are shown.
d
time between reheating events becomes shorter, and heat begins to build up.
te
Well away from either the internal stripe boundary or the free edge of the item, measurements from figure 14(a) indicate that the un-remelted width of the pool perpendicular to the scan vectors achieves a steady state value
Ac ce p
138.1 ± 10.0 µm, approximately equal to the expected nominal spacing between
the parallel scanning vectors of 140 µm, which is also the value that would be predicted by the Rosenthal model. However, near both the internal stripe boundary and the free surface, the melt pool width is larger for vectors which are moving away from the boundaries, sometimes by a significant amount. This locally increased width in turn partially covers the adjacent, previously processed vector making the apparent width of vectors approaching the boundary appear narrower than nominal, creating a pattern of alternating wide and narrow tracks. Figure 14(b) shows model output, with predicted thermal histories postprocessed to identify the un-remelted width of melt pools on the top surface at
35
Page 35 of 48
ip t cr us an M d te Ac ce p
Figure 14: Region of interest identified in figure 13. (a) BSE image showing spatial variation in the un-remelted apparent width of the scan vectors, along with scan vectors colored by time processed, and measurement location lines at 0.3, 1.5, and 3.0 mm [red, green, blue respectively] from the stripe boundary [black]. (b) DSM prediction of un-remelted top surface with the same vectors and measurement lines. Larger versions are available in the supplemental information.
36
Page 36 of 48
Experiment DSM
200 150 50
−100
Line 1: 0.3mm
−150 100
us
50 0 −50 −100
Line 2: 1.5mm
an
−150 50 0 −50 −100
Line 3: 3.0mm
−150 5
10
Scan Vector
15
20
te
d
0
M
Deviation from Ideal Unremelted Width [μm]
0 −50
cr
ip t
100
Ac ce p
Figure 15: Measurements of the un-remelted top surface in the region shown in figure 14 as predicted by the DSM (open symbols) and measured from the experiment (closed circles), expressed in terms of deviation from the hatch spacing, 140 µm, which is the prediction of the Rosenthal solution. The grey band indicates typical variation in observed width well away from a stripe boundary or free surface of ±10.0 µm.
z = 0, and figure 15 shows the un-remelted width measurements from both the experimental image as well as the model output along measurement lines labeled 1,2,3. The scan vectors have been discretized in the same fashion as employed in all earlier examples. The width is expressed in terms of deviation from the nominal 140 µm hatch spacing for each scan vector crossed by the measurement
line, beginning with vector 1 nearest to the free surface, and counting into the part. For Line 1 [red] at 300 µm from the stripe boundary, both the observed and modeled widths differ significantly from the nominal 140 µm hatch spacing, systematically alternating between narrower and wider. As described above, 37
Page 37 of 48
vectors moving away from the stripe boundary experience a pre-heat effect due to the previously scanned vector. They partially, or in some cases completely,
ip t
cover the previously deposited adjacent track, leading to the alternating wide-
narrow-wide behavior. Measurements in figure 15 show that the DSM also
cr
captures this trend, though the amplitude of the oscillation is under-predicted. Additionally, the vectors closest to the free surface on Line 1 show that as
us
the corner is approached, the overall scan vector length decreases, decreasing the time between re-heating events, increasing preheat, and therefore width. Qualitatively, it can also be seen that both the experiment and DSM predict
an
temporal overlap of melt pools from adjacent vectors, for the vectors closest to the intersection of the stripe-boundary and free surface. Further from the
as the beam is turning around.
M
corner, the melt pool induced by a particular track is able to completely solidify
Similar trends are observed for measurements along Line 2 [green] at 1.0 mm from the stripe boundary, though the amplitude of the alternating width is
d
weaker, as the pre-heat variation is generally smaller, again, with the DSM
te
under predicting the variation. Finally, for Line 3 [blue] at 3.0 mm from the stripe boundary, the observed and predicted width are both generally close to
Ac ce p
the expected value of 140 µm, with deviations apparent only for scan vectors close to the free-surface. Variations in the width observed along Line 3 in the DSM result are generally smaller than the experimentally observed variation in track width of ±10 µm in regions far from the stripe boundary or free surface. While the DSM is able to capture a number of trends in the spatial variation
of melt pool shape and duration molten, there are still discrepancies with the empirical results. Fluid flow, particularly near the beginning of a scan vector, is likely important in determining the extent to which the pool flows over the top of the adjacent track, with the expectation that flow driven by recoil pressure would drive liquid to recover adjacent tracks to a higher degree than predicted by the DSM. Latent heat of fusion was observed to have a slight impact on single track widths, but a more significant effect on pool lengh, and therefore duration molten. Longer lived molten pools would have more of an opportu38
Page 38 of 48
nity to overlap with those from adjacent tracks, again leading to a larger molten pool than predicted by the DSM. Additionally, the actual boundary condition is
ip t
somewhere between zero flux and uniform conductivity. However, similar calculations without the no-flux condition do not show a noted change in unremelted
cr
width for the end of vector features along the stripe-boundary, and thus is likely not the primary cause of discrepancy. Finally, we employ a uniform value of
us
T0 = 100 ◦C. We expect in practice that this value will have some spatial variation, however the inter-layer times are significant and the geometry of this build
an
simple enough that it is unlikely to drive the observed discrepancies. 6. Performance and Scaling
The DSM has been implemented as a filter within the DREAM.3D software
M
environment[42]. Specifically, the model has been implemented in C++ with thread parallelization performed using Intel’s Threaded Building Blocks (TBB). In an effort to quantify the practical scalability of the model, we performed
d
both a strong and weak scaling study. The study was carried out on an HP
te
Z840 with dual Intel Xeon E5-2687W V4 CPUs with base frequency 3.0 GHz, resulting in 48 logical processors. The strong scaling study involved evaluating
Ac ce p
3.75M spatial locations utilizing nLP =1, 6, 12, 24 and 48 logical processors. The weak scaling study utilized the same nLP values, but adjusted the number
of spatial evaluation locations Nj between 0.3125M, 1.875M, 3.75M, 7.5M and 15M correspondingly. The results of the two studies are listed in table 5, both in actual wall time and in scaling efficiency. For the weak scaling study, perfect load balancing was not explicitly enforced, rather TBB’s default load balancing scheme was used. The specific problem geometry and model parameters used are the same as described in section 5. The time resolution of the calculated thermal histories at each spatial evaluation location was selected to be 100 µs in all cases. As described in Section 2.5, the source discretization and output time resolution are held constant, and thus each spatial output location has, on average, the same number of evaluations of equation 9.
39
Page 39 of 48
Both strong and weak scaling efficiencies indicate that the DSM scales well. In both cases, scaling efficiency of ≈40% is achieved for the extreme cases. Fur-
ip t
thermore, while untested here, efficient scaling of parallelization in distributed
memory environments is also expected due to the uncoupled nature of the model,
cr
and is a means to further increase feasible values of Nj .
The case with Nj = 15 × 106 and nLP = 48 evaluated approximately 6300
us
output points per second of wall time on the hardware employed. Accounting for a 3 dimensional distribution of source points at the same density employed for the calculation in figure 14, using Rcut = 1 mm, the DSM could screen
an
temperature histories at points separated by 1 mm filling a cube with edge length 25 mm in approximately 60 s of wall time. Alternately, filling a cube 250 mm on edge (approximately a full build volume) with the same 1 mm grid of points
M
would require 15 h. For comparison, we estimate that an actual build of this height would require approximately 140 h of printing time using a generous estimate of 60 s per layer. Additionally, we note that generally significantly
d
less than 50% of build volumes are typically occupied with parts, and thus the
above.
te
spatial resolution could be increased while maintaining the same wall time cost
Ac ce p
There are opportunities to sample the scan vectors and place source points more efficiently than the uniformly spaced approach taken here with little or no loss of accuracy. Future work will focus on characterizing and quantifying such schemes, but initial estimates indicate a further reduction in the above mentioned 15 h calculation to approximately 2 h. Finally, while the above examples place the evaluation points on a cubic grid for simple comparison, we note again that there is complete freedom in placement and number of these points within the volume of interest, with no loss of accuracy for any particular thermal history. 7. Conclusion A discrete source model for powder bed fusion thermal history has been developed and verified against analytic calculations for single track geometries. 40
Page 40 of 48
6 twall E 1121 94 2167 95 -
12 twall 1179 -
E 89, 87 -
24 twall 971 1712 -
E 53 62 -
48 twall 655 2365
E 39 45
ip t
1 twall E 1053 100 12378 100 -
cr
nLP → N j × 10−6 0.3125 1.875 3.75 7.5 15
us
Table 5: Wall Time (twall ) and scaling efficiency (E) results for a strong and weak scaling study using several numbers of logical processors nLP . Wall Times are in seconds, efficiencies are in percents, and numbers of spatial evaluation points Nj are in millions. Note the two efficiencies listed for the nLP = 12 case are for the weak and strong scaling, respectively.
an
Convergence studies have been presented, and guidelines for choosing parameters to ensure required levels of accuracy have been developed, and an understanding of the range of approximation errors have been quantified.
M
A model calibration methodology using the geometry of experimental singletrack deposits has been presented, and it was shown that the DSM can be calibrated to match both the melt pool width and depth for Ti-6Al-4V.
d
The DSM was utilized to predict several spatial and temporal transient ef-
te
fects due to thermal interactions induced by multi-track deposits. Scan pattern and component geometry combine to produce spatial variation in pre-heat, and
Ac ce p
thus track geometry. Several reasons for quantitative discrepancies between experiment and DSM predictions have been presented. The DSM is an effective means of quickly generating an estimate of the local
thermal history induced by complex scan strategies when combined with arbitrary component geometry. We expect this tool will be useful for screening and identification of situations leading to anomalous thermal history, and ultimately to direct usage of higher fidelity models and experimental capabilities. Simple temporal discretization schemes for the energy source have been utilized here for clarity, but more optimized sampling for generation of source points in part scale calculations is expected to further enhance performance and is currently being pursued. Additionally, the thermal histories resulting from the present model are still cumbersome to work with, and reduced dimensional representations of
41
Page 41 of 48
thermal history will be a thrust of future investigation.
ip t
Acknowledgments
We wish to acknowledge Jonathan Miller at the Air Force Research Labo-
cr
ratory for many helpful conversations. Thomas Spears of GE Aviation and Joy
Gockel of Wright State University were instrumental in executing single track
us
deposition experiments, and Sonya Boone of UES collected images of multi-track deposits.
an
References
[1] T. J. Horn, O. L. Harrysson, Overview of current additive manufacturing technologies and selected applications, Science progress 95 (3) (2012) 255–
M
282.
[2] W. E. Frazier, Metal Additive Manufacturing: A Review, Journal of
d
Materials Engineering and Performance 23 (6) (2014) 1917–1928. doi:
te
{10.1007/s11665-014-0958-z}.
[3] Y. Huang, M. C. Leu, J. Mazumder, A. Donmez, Additive Manufacturing:
Ac ce p
Current State, Future Potential, Gaps and Needs, and Recommendations, Journal of Manufacturing Science and Engineering - Transactions of the ASME 137 (1). doi:{10.1115/1.4028725}.
[4] W. E. King, A. T. Anderson, R. M. Ferencz, N. E. Hodge, C. Kamath, S. A. Khairallah, A. M. Rubenchik, Laser powder bed fusion additive manufacturing of metals; physics, computational, and materials challenges, Applied Physics Reviews 2 (4). doi:{10.1063/1.4937809}.
[5] R. Jiang, R. Kleer, F. T. Piller, Predicting the future of additive manufacturing: A delphi study on economic and societal implications of 3d printing for 2030, Technological Forecasting and Social Change 117 (Supplement C) (2017) 84 – 97. doi:10.1016/j.techfore.2017.01.006.
42
Page 42 of 48
[6] W. King, A. T. Anderson, R. M. Ferencz, N. E. Hodge, C. Kamath, S. A. Khairallah, Overview of modelling and simulation of metal powder bed
ip t
fusion process at Lawrence Livermore National Laboratory, Materials Science and Technology 31 (8) (2015) 957–968. doi:10.1179/1743284714Y.
cr
0000000728.
[7] M. J. Matthews, G. Guss, S. A. Khairallah, A. M. Rubenchik, P. J. Depond,
us
W. E. King, Denudation of metal powder layers in laser powder bed fusion processes, Acta Materialia 114 (2016) 33–42. doi:{10.1016/j.actamat.
an
2016.05.017}.
[8] S. Ly, A. M. Rubenchik, S. A. Khairallah, G. Guss, M. J. Matthews, Metal vapor micro-jet controls material redistribution in laser powder bed fusion
s41598-017-04237-z.
M
additive manufacturing, Scientific Reports 7 (2017) 4085. doi:10.1038/
[9] S. A. Khairallah, A. Anderson, Mesoscopic simulation model of selec-
d
tive laser melting of stainless steel powder, Journal of Materials Process-
2014.06.001.
te
ing Technology 214 (11) (2014) 2627–2636. doi:10.1016/j.jmatprotec.
Ac ce p
[10] S. A. Khairallah, A. Anderson, A. M. Rubenchik, J. Florando, S. Wu, H. Lowdermilk, Simulation of the main physical processes in remote laser penetration with large laser spot size, AIP Advances 5 (4) (2015) 047120. doi:10.1063/1.4918284.
[11] S. A. Khairallah, A. T. Anderson, A. Rubenchik, W. E. King, Laser powderbed fusion additive manufacturing: Physics of complex melt flow and formation mechanisms of pores, spatter, and denudation zones, Acta Materialia 108 (2016) 36–45. doi:{10.1016/j.actamat.2016.02.014}.
[12] B. A. Cowles, D. G. Backman, R. E. Dutton, Update to recommended best practice for verification and validation of ICME methods and models for aerospace applications, Integrating Materials and Manufacturing Innovation 4 (1). doi:10.1186/s40192-014-0030-8. 43
Page 43 of 48
[13] D. Rosenthal, The theory of moving sources of heat and its application to
ip t
metal treatments, ASME, 1946. [14] M. Tang, P. C. Pistorius, J. L. Beuth, Prediction of lack-of-fusion porosity
39 – 48. doi:10.1016/j.addma.2016.12.001.
cr
for powder bed fusion, Additive Manufacturing 14 (Supplement C) (2017)
[15] J. Raplee, A. Plotkowski, M. M. Kirka, R. Dinwiddie, A. Okello, R. R.
us
Dehoff, S. S. Babu, Thermographic Microstructure Monitoring in Electron Beam Additive Manufacturing, Scientific Reports 7 (2017) 43554. doi:
an
10.1038/srep43554.
[16] S. P. Narra, R. Cunningham, J. Beuth, A. D. Rollett, Location specific solidification microstructure control in electron beam melting of ti-6al-4v, Ad-
M
ditive Manufacturing 19 (2018) 160–166. doi:{10.1016/j.addma.2017. 10.003}.
d
[17] E. R. Denlinger, J. C. Heigel, P. Michaleris, T. A. Palmer, Effect of interlayer dwell time on distortion and residual stress in additive manufacturing
te
of titanium and nickel alloys, Journal of Materials Processing Technology
Ac ce p
215 (2015) 123–131. doi:{10.1016/j.jmatprotec.2014.07.030}. [18] T. Eagar, N. Tsai, Temperature-Fields Produced by Traveling Distributed Heat-Sources, Welding Journal 62 (12) (1983) S346–S355.
[19] J. Goldak, A. Chakravarti, M. Bibby, A new finite element model for welding heat sources, Metallurgical Transactions B 15 (2) (1984) 299–305. doi:10.1007/BF02667333.
[20] N. Nguyen, A. Ohta, K. Matsuoka, N. Suzuki, Y. Maeda, Analytical solutions for transient temperature of semi-infinite body subjected to 3-D moving heat sources, Welding Journal 78 (8) (1999) 265S–274S. [21] N. Nguyen, Y. Mai, S. Simpson, A. Ohta, Analytical approximate solution for double ellipsoidal heat source in finite thick plate, Welding Journal 83 (3) (2004) 82S–93S. 44
Page 44 of 48
[22] M. Van Elsen, M. Baelmans, P. Mercelis, J. P. Kruth, Solutions for modelling moving heat sources in a semi-infinite medium and applications to
ip t
laser material processing, International Journal of Heat and Mass Trans-
fer 50 (23-24) (2007) 4872–4882. doi:{10.1016/j.ijheatmasstransfer.
cr
2007.02.044}.
[23] M. Akbari, D. Sinton, M. Bahrami, Moving heat sources in a half space: Ef-
us
fect of source geometry, Proceedings of the ASME 2009 Heat Transfer Summer Conference HT2009 (2009) 685–694doi:{10.1115/HT2009-88562}.
an
[24] R. Komanduri, Z. B. Hou, Thermal analysis of the arc welding process: Part i. general solutions, Metallurgical and Materials Transactions B 31 (6) (2000) 1353–1370. doi:10.1007/s11663-000-0022-2.
M
[25] R. Komanduri, Z. Hou, Thermal analysis of the laser surface transformation hardening process, International Journal of Heat and Mass Transfer 44 (15)
d
(2001) 2845–2862. doi:10.1016/S0017-9310(00)00316-1. [26] R. Komanduri, Z. B. Hou, Thermal analysis of the arc welding process:
te
Part ii. effect of variation of thermophysical properties with temperature, Metallurgical and Materials Transactions B 32 (3) (2001) 483–499. doi:
Ac ce p
10.1007/s11663-001-0034-6.
[27] J. Winczek, Analytical solution to transient temperature field in a halfinfinite body caused by moving volumetric heat source, International Journal of Heat and Mass Transfer 53 (25-26) (2010) 5774–5781. doi: {10.1016/j.ijheatmasstransfer.2010.07.065}.
[28] S. Mohanty, J. Hattel, Cellular scanning strategy for selective laser melting Capturing thermal trends with a low-fidelity, pseudo-analytical model, Mathematical Problems in Engineering 2014 (2014) 715058. doi:10.1155/ 2014/715058. [29] S. Mohanty, J. H. Hattel, Cellular scanning strategy for selective laser melting: Generating reliable, optimized scanning paths and processing param45
Page 45 of 48
eters, Proc. SPIE Volume 9353, Laser 3D Manufacturing II 9353 (2015)
ip t
93530U. doi:10.1117/12.2079957. [30] S. Mohanty, J. H. Hattel, Reducing residual stresses and deformations in
selective laser melting through multi-level multi-scale optimization of cel-
III 9738 (2016) 93780Z. doi:10.1117/12.2212490.
cr
lular scanning strategy, Proc. SPIE Volume 9738, Laser 3D Manufacturing
us
[31] A. Plotkowski, M. Kirka, S. Babu, Verification and validation of a rapid heat transfer calculation methodology for transient melt pool solidification
an
conditions in powder bed metal additive manufacturing, Additive Manufacturing 18 (Supplement C) (2017) 256 – 268. doi:10.1016/j.addma. 2017.10.017.
M
[32] C. D. Boley, S. A. Khairallah, A. M. Rubenchik, Calculation of laser absorption by metal powders in additive manufacturing, Applied Optics 54 (9)
d
(2015) 2477–2482. doi:10.1364/AO.54.002477. [33] A. Rubenchik, S. Wu, S. Mitchell, I. Golosker, M. Leblanc, N. Peterson,
te
Direct measurements of temperature-dependent laser absorptivity of metal powders, Applied Optics 54 (24) (2015) 7230–7233. doi:10.1364/AO.54.
Ac ce p
007230.
[34] J. Trapp, A. M. Rubenchik, G. Guss, M. J. Matthews, In situ absorptivity measurements of metallic powders during laser powder-bed fusion additive manufacturing, Applied Materials Today 9 (Supplement C) (2017) 341 – 349. doi:10.1016/j.apmt.2017.08.006.
[35] S. Lu, H. Fujii, K. Nogi, Sensitivity of Marangoni convection and weld shape variations to welding parameters in O-2-Ar shielded GTA welding, Scripta Materialia 51 (3) (2004) 271–277. doi:10.1016/j.scriptamat. 2004.03.004.
[36] K. C. Mills, Recommended Properties of Commercial Alloys, Woodhead Publishing Ltd. & ASM International, 2002. 46
Page 46 of 48
[37] R. R. Dehoff, M. M. Kirka, W. J. Sames, H. Bilheux, A. S. Tremsin, L. E. Lowe, S. S. Babu, Site specific control of crystallographic grain orientation
ip t
through electron beam additive manufacturing, Materials Science and Technology 31 (8) (2015) 931–938. doi:10.1179/1743284714Y.0000000734.
cr
[38] R. R. Dehoff, M. M. Kirka, F. A. List, III, K. A. Unocic, W. J. Sames, Crystallographic texture engineering through novel melt strategies via electron
us
beam melting: Inconel 718, Materials Science and Technology 31 (8) (2015) 939–944. doi:10.1179/1743284714Y.0000000697.
an
[39] D. G. Duffy, Green’s Functions with Applications, Second Edition, Taylor & Francis Group, 2015.
[40] J. D. Hunt, STEADY-STATE COLUMNAR AND EQUIAXED GROWTH
M
OF DENDRITES AND EUTECTIC, Materials Science and Engineering 65 (1) (1984) 75–83. doi:10.1016/0025-5416(84)90201-5.
d
[41] W. Kurz, C. Bezencon, M. Gaumann, Columnar to equiaxed transition
te
in solidification processing, Science and Technology of Advanced Materials 2 (1, SI) (2001) 185–191. doi:10.1016/S1468-6996(01)00047-X.
Ac ce p
[42] M. A. Groeber, M. A. Jackson, Dream.3d: A digital representation environment for the analysis of microstructure in 3d, Integrating Materials and Manufacturing Innovation 3 (1) (2014) 5. doi:10.1186/2193-9772-3-5.
[43] R. Balluffi, S. Allen, W. Carter, Kinetics of Materials, John Wiley & Sons, 2005.
[44] H. Hu, S. A. Argyropoulos, Mathematical modelling of solidification and melting: a review, Modelling and Simulation in Materials Science and Engineering 4 (4) (1996) 371. [45] V. R. Dave, D. A. Hartman, M. J. Cola, In-Process Quality Assurance for Aerospace Welding, Welding Journal 88 (2) (2009) 28–32.
47
Page 47 of 48
[46] M. Uchic, M. Groeber, M. Shah, P. Callahan, A. Shiveley, M. Scott, M. Chapman, J. Spowart, An Automated Multi-Modal Serial Sectioning
ip t
System for Characterization of Grain-Scale Microstructures in Engineering
Materials, in: DeGraef, M and Poulsen, HF and Lewis, A and Simmons, J
cr
and Spanos, G (Ed.), Proceedings of the 1st International Conference on 3D Materials Science, 2012, pp. 195–202, Seven Springs, PA, Jul. 08-12,
Ac ce p
te
d
M
an
us
2012.
48
Page 48 of 48
S2214-8604(18)30778-4 https://doi.org/doi:10.1016/j.addma.2018.12.004 ADDMA 603
To appear in: Received date: Revised date: Accepted date:
1 October 2018 7 November 2018 3 December 2018
Please cite this article as: Edwin J. Schwalbach, Sean P. Donegan, Michael G. Chapman, Kevin J. Chaput, Michael A. Groeber, A Discrete Source Model of Powder Bed Fusion Additive Manufacturing Thermal History, (2018), https://doi.org/10.1016/j.addma.2018.12.004 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ip t
A Discrete Source Model of Powder Bed Fusion Additive Manufacturing Thermal History Edwin J. Schwalbacha,∗, Sean P. Donegana , Michael G. Chapmanb , Kevin J. Chaputa , Michael A. Groebera,c,∗ a Air
an
us
cr
Force Research Laboratory, Materials and Manufacturing Directorate, Wright-Patterson Air Force Base, OH 45433 b UES Inc. Beavercreek, OH 45432 c Ohio State University, Department of Integrated Systems Engineering, Columbus, OH 43210
Abstract
Significant attention has been focused on modeling of metallic additive manu-
M
facturing (AM) processes, with the initial aim of predicting local thermal history, and ultimately structure and properties. Existing models range greatly in physical complexity and computational cost, and the implications of various
d
simplifying assumption often go unassessed. In the present work, we first formu-
te
late a fast acting Discrete Source Model (DSM) capable of handling the complex processing often encountered in metal powder bed fusion AM. We then assess
Ac ce p
implications of the source representation, details of the numeric implementation, as well as effects of boundary conditions and thermophysical parameters. We verify the DSM implementation against simple numerical thermal predictions, calibrate it with single track deposit experiments, validate outputs against multitrack deposits, and finally quantify the scaling performance. The DSM is an effective means of quickly generating an estimate of the local thermal history induced by complex scan strategies when combined with arbitrary component geometry. While a number of approximations limit its quantitative accuracy, the inexpensive nature and ability to treat complex processing plans suggests it will be useful for screening and identification of regions experiencing anomalous ∗ Corresponding
author Email addresses: [email protected] (Edwin J. Schwalbach), [email protected] (Michael A. Groeber)
Preprint submitted to Elsevier
November 7, 2018
Page 1 of 48
thermal history. Such a capability is necessary to direct usage of higher fidelity, more expensive models and experimental resources.
ip t
Keywords: selective laser melting, powder bed fusion, thermal modeling,
cr
model verification and validation
1. Introduction
us
Additive Manufacturing (AM) technologies have received considerable research and development attention in the last decade[1, 2, 3, 4, 5]. This collection of materials processing technologies has been touted to provide designers
an
with a number of potential benefits, including manufacturability of unique geometries offering significant improvement in component performance, as well
M
as reductions in either manufacturing cost, lead time, or both. The concept of processing-structure-property relationships is not unique to materials fabricated via AM means, but AM processes can have considerable complexity and offer an
d
extreme degree of location specificity compared to many conventional processes. Therefore, much attention has been focused on modeling of AM processes, with
te
the aim of initially predicting local structure and ultimately properties. With respect to metals AM processes, the local thermal history induced by
Ac ce p
spatially and temporally complex energy input procedures affects microstructure, defects, residual stresses, and ultimately material performance. In this work, we focus on predicting the thermal histories generated by powder bed fusion (PBF) AM processes. Such models have been the focus of intense study over the past several years[6]. Depending on the specific problem, the range of physical phenomena and energy transport mechanisms involved can be extensive: conduction, radiation, phase transformations, and convective transport are all active, with the later driven by capillarity, Marangoni forces, bouyancy, and even recoil pressure[7, 8]. Models incorporating various combinations of such effects have generated considerable insight into the detailed transport phenomena that drive thermal field evolution[9, 10, 11], though high-fidelity solutions to such multiphysics problems are currently still expensive to compute and hence
2
Page 2 of 48
limited in terms of the volume of material they can interrogate. Models that employ reasonable simplifying assumptions in order to efficiently
ip t
incorporate the relevant physics are desirable for screening and trend prediction
purposes, particularly in light of the the widening AM processing parameter
cr
space. Such models inherently rely on some degree of calibration, whether to empirical data or results from higher fidelity models, and must be used with
us
care. It is critical to understand the application intent, range of applicability, and supporting data requirements for such models [12].
The so-called Rosenthal solution is a model originally developed to analyze
an
fusion welding[13], though it has seen much attention from the metal AM community. This is a computationally inexpensive analytic representation of the temperature field induced by a point source of energy moving in a straight line
M
at a constant velocity across a semi-infinite body. Furthermore, it accounts for only conductive transport, and assumes temperature independent thermophysical parameters. As can be inferred from this list of assumptions, this model
d
has several limitations. First, as the energy input is a point source, there is a
te
singularity in the resulting temperature field at the source location. The substrate is taken to be semi-infinite, which may not be unreasonable for very short
Ac ce p
time scales, though it is an increasingly poor assumption as elapsed time becomes long compared the to the thermal diffusion time scale associated with component size. Additionally, the assumption that thermal conduction is the only relevant transport mechanism is questionable in the immediate vicinity of the source and resultant melt pool, where we have already noted that radiation, convective transport, and evaporation can be significant. Finally, and perhaps most importantly, because the energy source moves uniformly in space and time, the model cannot directly describe transient effects induced by source direction changes or modulations in processing conditions such as power or speed, or interactions of multiple sources. In spite of these limitations, this model is useful when operating in the socalled conduction mode. It predicts empirically observed trends in melt pool size and aspect ratio and can be calibrated with an efficiency factor decrementing 3
Page 3 of 48
the input power to bring model output into reasonable quantitative agreement with empirical experience over limited operating parameter ranges. Tang, Pis-
ip t
torius and Beuth derive a melt pool width from the canonical Rosenthal thermal
(1)
us
top surface of the deposit is approximately 1/2 8P w= , eπρcp (Tm − T0 ) v
cr
field[14]. Equation 20 from that work indicates the melt pool width w on the
where P is the absorbed power, ρ is the density, cp is the mass specific heat capacity, Tm is the melting temperature, T0 is the spatially uniform starting or
an
initial temperature of the substrate, and v is the velocity of the beam. Similarly, the length measured from the instantaneous center of the energy source to the trailing edge of the melt pool is found to be
P , 2πκ (Tm − T0 )
M
l=
(2)
where κ is the thermal conductivity. For the case of single track deposits,
te
characterized.
d
all of the parameters appearing in equations 1 and 2 can be reasonably well
In cases where more complex energy input strategies are employed, as is com-
Ac ce p
mon in PBF, the intial temperature T0 can vary significantly with position[15], which in turn has been observed to affect microstructure[16] and residual stress [17]. In principle, the Rosenthal solution could be utilized to estimate thermal field details for cases with more complex energy input strategies, assuming the starting temperature T0 could be easily predicted, and with an additional assumption that the temperature in the local vicinity of the point of interest is sufficiently uniform. However, this is particularly challenging as the energy input path becomes complex, as in the case for typical PBF scan strategies. A number of other researchers have expanded on the model of Rosenthal by replacing the point source with various spatially distributed energy sources [18, 19, 20, 21, 22, 23], effectively regularizing the temperature discontinuity. Additionally, others have focused on expanding the various limitations on path complexity [24, 25, 26, 27, 28, 29, 30, 14, 31]. The main objective in the fol4
Page 4 of 48
lowing work is to not only develop a computationally inexpensive model to predict transient thermal fields produced by arbitrary energy input procedures
ip t
relevant to current and future PBF AM processes, but furthermore to assess model accuracy, establish ranges of applicability, and quantify computational
cr
performance.
The remainder of the paper is organized as follows: We discuss relevant phys-
us
ical assumptions and limitations, formulate a transient thermal history model, and then outline an efficient numeric implementation. We then demonstrate model convergence and validate model output against analytic solutions un-
an
der simple conditions. We calibrate model output against simple experimental results, and then validate model predictions against experimental observations from more complex multipass geometric configurations. Finally, we con-
M
clude with a discussion of model performance and scaling, and then outline future model use in conjunction with reduced order thermal history representations. Throughout this work, we will employ thermophsyical parameters for
d
the aerospace relevant alloy Titanium-6wt%Aluminum-4wt%Vanadium (Ti-6Al-
te
4V). 2. Model Formulation
Ac ce p
2.1. Basic Assumptions
The objective of the present model is to predict an approximate temperature
field during PBF AM processes resulting from (potentially multiple) sources of energy input, such as a laser or electron beam (e-beam), moving along spatially complex paths over a component being fabricated. The primary mechanisms driving thermal field evolution are the energy input itself and resulting thermal conduction.
Energy fluxes due to radiation and latent heat of evaporation can also be significant. Their effects are strongest at very high temperatures and thus their primary influence will be approximately coincident with the molten pool, which is in turn often coincident with the energy input source in both space and time. Laser and e-beam sources can also have absorption efficiencies significantly less 5
Page 5 of 48
than 1.0 [32, 33], and are known to vary appreciably with processing conditions, in particular near the conduction to keyholing mode transition[34]. The present
ip t
model employs an efficiency factor η to decrement the nominal input power P
to account for both the absorption efficiency, as well as losses in the immediate
cr
vicinity of the source due to radiation and evaporation.
Convective flux within the liquid phase is also known to have a significant
us
effect on the shape of melt pool in certain regimes[35, 11]. Rather than directly accounting for these complex physics, the present model modifies the geometry of a volumetric energy source to reproduce the observed melt pool geometries,
an
in particular its depth to width ratio, similar to the approach of Goldak et al.[19]. This shape factor must also be determined via calibration. A process for calibration against experimental data will be discussed in detail in section 4.2,
M
however for both the efficiency and shape parameters, calibration could also be determined by comparison to trusted, higher fidelity models. The latent heat stored or released during liquid-solid and solid-state trans-
d
formations can also affect the thermal field. The net energy contributions from
te
such terms sum to zero over long enough time periods since the amounts of material melted and solidified by the AM PBF process are equivalent. How-
Ac ce p
ever, in the vicinity of the phase transformations these terms will modify the temperature distribution, particularly near the solid-liquid interface. The Stefan number measures the ratio of sensible to latent heat 1 St = Esens /Hsl = Hsl
ZTm cp (T ) dT,
(3)
T0
where Hsl is the mass specific enthalpy of the solid-liquid transformation and cp is the mass specific heat capacity. Heating Ti-6Al-4V from room temperature to melting, we find St = 3.81, using values from [36]. Thus, the sensible heat is roughly four times the latent heat of the transformation, and we expect latent heat to be a secondary, though not completely negligible, effect. We can further compare this to the total amount of energy the beam delivers to the melt pool. Approximating the melt pool volume V as a quarter sphere
6
Page 6 of 48
with radius w/2 and a quarter ellipsoid with axes a = b = w/2, and c = l, where
ip t
w and l are determined from equations 1 and 2, the volume of the melt pool is πw2 w V = +l . (4) 12 2
stored in the liquid is
(5)
us
Esens
ZTm ρ (T ) cp (T ) dT, =V
cr
Assuming the whole volume of this melt pool is at Tm , the total sensible heat
T0
(ignoring the latent heat of solid-state transformations which are typically smaller),
an
and the total latent heat of fusion stored in the melt pool Ef = V ρHsl . Finally, the total energy delivered to the melt pool by the beam is approximately ηP l . v
(6)
M
Eabs =
Substituting typical values for Ti-6Al-4V and assuming an efficiency η = 50% based on reports in references [32, 33, 34], we find Esens /Eabs = 0.275, and
d
Ef /Eabs = 0.071. The total energy delivered to the melt pool is therefore large compared to the energy required to bring it to the melting point, which is in turn
te
large compared to Ef . We will assess the implications of these circumstances further in Section 3.4.
Ac ce p
In light of the above assumptions and simplifications, energy balance in
the presence of N generic volumetric energy sources sˆi which are functions of position ~r and time t, the relevant energy transport equation is N
X ∂ (ρcp T ) = ∇κ · ∇T + sˆi (~r, t) , ∂t i
(7)
where T is temperature, ρ is mass density, cp is mass specific heat capacity, κ is thermal conductivity, and sˆi is power per unit volume introduced by the ith source. With the further assumption that the thermal and transport properties are independent of temperature and uniform in space and time, the above reduces to
N
X sˆi (~r, t) ∂T = α∇2 T + , ∂t ρcp i
(8)
where the thermal diffusivity α is κ/ (ρcp ). 7
Page 7 of 48
2.2. Energy Source Representation
ip t
For conventional PBF AM processes there is a single continuously scanning energy source. To simplify the expression of such a problem, particularly for
cases where the beam moves in arbitrary spatial patterns, we instead represent
cr
this moving source as a series of discrete, stationary energy source terms. The
details of sˆi (~r, t) in equation 8 are then determined by sampling the position of
us
the continuously moving beam, for now at a fixed rate ∆t, and placing stationary sources at these locations ~ri = hxi , yi , zi i, and then activating them at the time τi at which the moving beam passes over them. We refer to this as the Discrete
an
Source Model (DSM), and will demonstrate that on time scales greater than ∆t the thermal profile produced by such a model will converge to that of the original continuously moving source. An example single-track vector is shown
interrogation point of interest.
M
in figure 1(a), with its DSM representation in figure 1(b) near a particular
While formulated in the context of a scanning beam type energy input mode,
d
with appropriate choices of source positions and activation times this method
te
can be used to represent an arbitrary energy input process. It is not limited to a single beam, and in fact could also account for shaped sources such as lines
Ac ce p
or areas. In some cases such as pulsed lasers or the “spot melting” approach employed in [37, 38] to control texture and grain morphology, the discrete approach described here is actually a faithful representation of the energy input process.
The individual source terms employed in the DSM are denoted si and have
a Gaussian distribution in space centered at locations ~ri and supply an impulse of energy at time τi . The effect of a source at location ~ri on the temperature at arbitrary position ~rj and time t is si (Rij , t) =
Ai δ (t − τi ) 3/2
(2πσ 2 )
2 Rij exp − 2 2σ
! .
(9)
Here, Ai δ (t − τi ) is the power of the source with δ the Dirac delta function, σ sets the physical size over which this energy is deposited, and the distance
8
Page 8 of 48
ip t cr 5
10
15
20
25
X[mm]
30
0.025 0.000 −0.025 19.7
19.8
0.0152
(c)
19.9
0.0153
7 5
20.2
20.3
0.0155
0.0156
4
te
Y[mm]
20.1
0.0154
τi[s]
40
d
6
20.0
X[mm]
35
an
(b)
Y[mm]
0
us
1 0 −1
M
Y[mm]
(a)
3 2 1
Ac ce p
0
0
2
4
X[mm]
6
Figure 1: Schematic geometries considered for various problems. (a) 40 mm single vector showing original scan vector (grey arrow) prior to discretization and the interrogation point (black circle). (b) Zoom of the single bead in the region of the interrogation (black circle) showing discrete source locations r~i (colored circles) colored by their timings τi . Open black circles indicate 4σ diameters of two specific sources. (c) 7 × 7 mm square with alternating vectors and several interrogation points.
9
Page 9 of 48
Rij = |~rj − ~ri | is defined for convenience. Note that the units of Ai δ (t) are energy per unit time, and thus the units of si are energy per unit time per
ip t
unit volume. In this model, the source deposits all of it’s energy at t = τi , and
then becomes inactive, though the energy deposited continues to influence the
cr
temperature distribution for all t ≥ τi .
Integrating each si (Rij , t) over the full domain in space and time gives the
us
total amount of energy contributed by the ith source, and equating this total to the effect of the continuous source with efficiency η, τi + ∆t 2
Z Z
Z
ηP (t) dt ≈ ηi Pi ∆t
(10)
an
si (Rij , t) d~rj dt = t
V
τi − ∆t 2
where we have assumed that the continuous source varies in power slowly com-
M
pared to the time scale ∆t and that the effective power of the source ηi Pi is simply equal to the actual power evaluated at τi , ηP (τi ). This approximation is exact when the power is constant over τi − ∆t/2 < t < τi + ∆t/2.
d
For a source located on the object’s top surface defined by zi = 0, we inte-
te
grate over the half space V : −∞ < x < ∞, −∞ < y < ∞, z ≤ 0 and over all time to find the total energy deposited by source i:
Ac ce p
ηi P i
Ai (2πσ 2 )3/2
∆t = R R2 δ (t − τi ) V exp − 2σij2 d~rj dt, −∞
(11)
R∞
Executing the integral and rearranging, we have Ai = 2ηi Pi ∆t.
(12)
The factor of two comes from the fact that half of the total energy introduced by the source is deposited above the plane z = 0. The strength is effectively doubled to ensure that the energy distributed below this plane of symmetry is equivalent to the true continuous source. This also effectively results in a zero-flux boundary condition along the top surface of the part, z = 0.
10
Page 10 of 48
Substituting equation 12 into the original transport equation, we have
+
PN i
= α∇2 T 2 Rij 2ηi Pi ∆t δ (t − τi ) . − 2 3/2 exp 2 2σ ρc (2πσ ) p
(13)
ip t
∂T ∂t
cr
As each new layer in a multilayer deposit is added the coordinate system is adjusted such that the instantaneous top surface of the deposit remains at
us
z = 0. In this reference frame, all sources activated during previous layers will have zi < 0. Such sources can directly influence only material that was added during or prior to the layer during which they became active. Combinations of
an
~ri , ~rj where zi < zj are excluded from the sum in equation 13. This effectively indicates that energy deposited by si cannot diffuse up towards positive z. For cases where the system cools significantly between layers this is not a limiting
M
assumption and the effect of all prior sources is easily captured by the uniform background temperature T0 . However, if the interlayer time is short, or if there is significant non-uniformity of temperature due to heating of newly deposited
d
powder by previously deposited energy, this assumption must be revised. One
te
method to do this would be to accommodate a non-uniform initial temperature condition determined as a function of sources from previous layers, though this
Ac ce p
would require marching the solution in time. 2.3. Thermal Solution
We utilize a Green’s function methodology to obtain solutions to equation
13 in an infinite medium with spatially uniform thermophysical properties and initial temperature T0 [39]. For the case of spherical Gaussian sources we find ! N 2 X Rij −3/2 T (~rj , t) = T0 + Ci Θ(t − τi )λi exp − , (14) 2λi i
where Θ is the Heaviside step function (Θ (x) = 0 if x < 0, and 1 otherwise), and for convenience we have defined Ci =
ηi Pi ∆t √ , ρ cp 2π 3/2
λi = σ 2 + 2α (t − τi ) .
(15) (16)
11
Page 11 of 48
An analogous approach can be taken for a source with an elliptical Gaussian energy distribution. With the ith source dimensions in the x, y, z directions
∆t , 2π 3/2 ρcp
λqi = σqi 2 + 2α (t − τi ) ,
(17)
q = x, y, z,
T (~rj , t) = T0 +
"
N X
(18)
us
the temperature is expressed as
cr
C=√
ip t
respectively being σxi , σyi , and σzi , and defining
Θ(t − τi ) √ ηi Pi C
λxi λyi λzi
i x 2 − 2λijxi
yij 2 2λyi
zij 2 2λzi
an
× exp
−
−
(19)
#
.
where qij = qj − qi for q = x, y, z. In the following work we employ the simplifi-
a plane normal to the beam.
d
2.4. Derivatives
M
cation σx = σy as most laser energy distributions are approximately circular in
~ and cooling rate ∂T /∂t = R are imporThe temperature gradient ∇T = G
te
tant conditions that influence solidification morphology [40, 41]. Solidification 2 ~ , with ~ G front velocity ~v is related to these quantities according to ~v = −RG/
Ac ce p
the sign convention such that the velocity vector is parallel to the thermal gradient in regions that are solidifying (R < 0), and anti-parallel on melting (R > 0). ~ is useful in identifying a position ~r which lies on Additionally, knowledge of G a particular isotherm of interest via Newton’s method. This technique is one method employed to find the position of points on the melting isotherm, and thus measure melt pool dimensions in the following. Analytic expressions for the temporal and spatial derivatives are directly
evaluated from equations 14 and 19, and are listed in the supplemental information for convenience. 2.5. Numeric Implementation Equation 14 indicates that the ith source affects the temperature at all points in space when t ≥ τi . However, the importance of the ith term at an arbitrary 12
Page 12 of 48
location ~rj depends on both the distance Rij and elapsed time t − τi since the source has been activated; increasing either distance or time results in dimin-
ip t
ished influence. Neglecting relatively weak terms in the sum will reduce the computational expense, and, if done with care, will have a minor and measur-
cr
able effect on solution accuracy.
Examining equations 14 and 19, it is clear that when Rij is sufficiently large
us
compared to the denominator of the exponential term, the source’s influence will be weak. Furthermore, σ 2 is typically small compared to the diffusive term as t − τi increases, suggesting that sources which satisfy the simplified condition
an
2 Rij 4α (t − τi ) are relatively unimportant. In practice, a finite cutoff radius
Rcut is chosen, and only terms satisfying the condition Rij ≤ Rcut are retained in the sums in equations 14 and 19. For common scan strategies source terms
M
are approximately uniformly distributed in space, and thus the number of terms retained in the sum will scale approximately as the volume contained within the cutoff distance ∝ Rcut 3 .
d
While the spatial extent of a source term’s influence increases with elapsed
te
time, the magnitude of its impact decreases. Setting the pre-exponential coefficient in equation 14 equal to a minimum temperature change of interest ∆Tmin
Ac ce p
and solving for elapsed time, we have t − τi = Ci / (2α ∆Tmin ). This relationship gives a rational methodology to select a time cutoff tcut , after which source i is no longer important; i is dropped from the sum if t−τi < tcut . We will assess the influence of these approximations on both execution time and solution accuracy in the following sections.
Note that equations 14 and 19 indicate that the computation of T (~rj , t) is
independent of the temperatures at any other points in space or time. Two major implications of this computational decoupling are: 1.) the computation of many space-time temperature outputs is embarrassingly parallel and, as will be seen, scales reasonably well across multiple processors, and 2.) the accuracy of the thermal field approximation is not dependent on the choice space-time output points (~rj , t), i.e. there is no grid stability or output mesh dependence of the solution. 13
Page 13 of 48
ip t cr
us
initialize Tjk = T0 identify active sources for evaluation points (~rj , t) for all evaluation points j [parallel] do for active sources i where Rij < Rcut do for times k where τi ≤ tk < τi + tcut do ∆Tijk = si (Rij , tk ) (equation 9) Tjk = Tjk + ∆Tijk end for end for end for
an
Figure 2: Pseudocode describing the parallel implementation of the present model. Parallelism occurs across the outermost for loop. Index i refers to discrete sources, j to the desired output spatial locations, and k to the desired time output points. Note that the choice of time output points does not have to be equivalent across all spatial output locations.
Parallelization is performed over the spatial evaluation locations, with each ~rj calculation looping only over the source locations ~ri satisfying the time and
M
space cutoffs tcut and Rcut respectively. Each source location evaluates equation 9 for all desired output times of interest. For the specific case of a constant
d
output time interval, the number of time output points is simply tcut divided by the desired output time resolution. Each evaluation of equation 9 is added
te
to a time point for the current spatial evaluation location. Pseudocode detailing this process is shown in figure 2. From the pseudocode,
Ac ce p
it is clear that reducing either the time resolution of evaluations or the source discretization will reduce the computational time of any given spatial evaluation iteration of the parallel loop. An implementation of this algorithm within the DREAM.3D software environment will be described in Section 6[42]. Specifically, a performance study with the source discretization and time resolution are held constant, resulting in each spatial evaluation location in each run of the model having, on average, the same number of evaluations of equation 9. 3. Model Verification In this section, we establish the impact of various physical and numerical approximations and verify the model by reproducing a known analytic result. Section 3.1 assesses the effect of the source discretization, and sections 3.2 and
14
Page 14 of 48
ip t
Description nominal power efficiency beam radius beam velocity melting temp. property temp. mass density heat capacity conductivity starting temp. cutoff distance source discretization cutoff time
cr
Unit W − µm mm s−1 ◦ C ◦ C kg m−3 J kg−1 K−1 W m−1 K−1 ◦ C mm µs s
us
Value 290 1.0 18.75 1300 1650 1203 4252 678.6 22.93 25.0 1.0 20 ∞
an
Symbol P η σ v Tm Tp ρ cp κ T0 Rcut ∆t tcut
M
Table 1: Default parameters utilized for all simulations, unless otherwise noted. The thermophysical parameters are representative of Ti-6Al-4V at a temperature of Tp = 1203 ◦C, as described in [36].
3.6 address numerical implementation effects. Section 3.4 deals with the implications of constant thermophysical parameters and neglect of latent heat.
d
Assumptions regarding boundary conditions and powder bed conductivity are
te
addressed in section 3.5. We defer effects of the neglect of fluid flow until discussing energy source shape in section 4.2.
Ac ce p
Table 1 indicates process, physical, and numeric parameters employed in calculations unless otherwise noted, and figure 1 shows relevant scan geometries that will be considered. Again, thermophysical parameters throughout are representative of Ti-6Al-4V. 3.1. Source Discretization
The relationship expressed in equation 10 guarantees that the same net
energy is deposited independent of the choice of ∆t. However, to accurately
represent the thermal history from a continuously scanning energy source, particularly in the temporal vicinity of an energy input event, an appropriate ∆t must be chosen. Excessively small ∆t will incur increased computational cost and large ∆t will induce spurious thermal excursions. A reasonable guideline to ensure accuracy in the vicinity of an energy input even is that ∆t must be
15
Page 15 of 48
sufficiently small with respect to the time tb = 4σ/v, (where σ = σx = σy in the case of the elliptical source) which is simply the elapsed time required for
ip t
the source to traverse a distance equal to its 4σ beam diameter. It is recognized
that this resolution is higher than needed to represent more spatially distant
cr
input events.
Figure 3 shows the temperature computed at a point on the center of the
us
40 mm long single-bead deposition vector shown in figure 1(a) and (b), as a function of time elapsed since the source was centered on the interrogation point at absolute time t0 . The calculations are performed using the conditions
an
indicated in Table 1, with the exception of the value of ∆t. The solutions exhibit temperature jumps as each discrete source is activated which are strongest in magnitude at times around t0 (see inset). These jumps decrease in magnitude
M
and the solutions become increasingly smooth as ∆t is decreased to a value well below tb . It is also clear that the solutions converge as the elapsed time becomes large compared to ∆t. For the conditions chosen in figure 3, all curves
d
are effectively equivalent approximately 1ms after the beam passes. Unless
te
otherwise stated, for the remainder of this work we employ ∆t = 2 × 10−5 s (or ∆t/tb ≈ 1/3). This discretization is shown in figure 1(b), along with black circles
Ac ce p
representing the spatial extent of two successively activated sources. While there are still visible thermal excursions in figure 3 for the chosen ∆t, the discrepancies
decay at temperatures that are still well above the melting point and thus will not have a significant impact on the predicted melt pool geometry. 3.2. Spatial Cutoff
Next we consider the effect of the spatial cutoff on the quality of the solution.
We employ the same single bead deposit geometry as described in the preceding section, but now vary the spatial cutoff distance Rcut . Figure 4 shows the temperature profile at the center of the bead as a function of time elapsed since the source was centered on the interrogation point at absolute time t0 . The
bounding black curve is the baseline solution with all terms retained in the sum (i.e. Rcut → ∞), and other curves show the result utilizing a range of finite Rcut .
16
Page 16 of 48
105
104 102 −50
0
50 t − t0[μs]
Δt/tb = 3.5 Δt/tb = 0.35 Δt/tb = 0.035 Δt/tb = 0.0035 100
150
103
250
500
750
1000 1250 t − t0[μs]
1500
1750
2000
an
0
ip t
Δt = 2e-04, Δt = 2e-05, Δt = 2e-06, Δt = 2e-07,
103
cr
104
us
T[C]
T[C]
105
d
M
Figure 3: Cooling curves for a point at the center of a single bead deposit using several values of the source discretization time scale ∆t. The beam passes over the point in question at t0 . Inset: expanded time scale surrounding the energy input event. Note: tb = 4σ/v = 58 µs is elapsed required for the beam to traversed it’s 4σ diameter. Nearly complete overlap of conditions 2 × 10−6 s and 2 × 10−7 s indicates convergence with respect to ∆t.
Dashed and dotted curves indicate temperature decay at rates proportional to
te
t−1 and t−3/2 respectively, as a guide to the eye. It is clear that the effect of the spatial cutoff is most siginicant for times
Ac ce p
well after the energy input event. The deviation between the baseline full solution and the various approximations are concentrated primarily at large elapsed times. The accuracy of the thermal history increases as the static cutoff distance is increased. With Rcut = ∞ the solution decays as t−1 for the duration of the calculation (10 s), which is characteristic of diffusive thermal decay in a 2D problem[43]. While the geometry in question is 3D in all cases, the motion of the beam along one direction effectively suppresses thermal diffusion in that dimension resulting in 2D cooling behavior. For cases with finite Rcut , as the elapsed time surpasses a value given by Rcut 2 / (4α) the absence of sources beyond Rcut is “felt” and the temperature decay transitions from cooling proportional to t−1
towards t−3/2 , characteristic of a 3D diffusion case. Ultimately, an appropriate value of Rcut must be chosen to balance compu-
17
Page 17 of 48
105
ip t
104
Rcut = ∞ Rcut = 0.1 [mm] Rcut = 0.3 [mm] Rcut = 1.0 [mm] Rcut = 3.0 [mm] Rcut = 10.0 [mm] ΔT ∝ t −1 ΔT ∝ t −32
101 100 10−1 10−2 10−6
10−5
10−4
10−3
10−2
t − t0 [s]
10−1
100
us
cr
102
101
an
T − T0 [C]
103
M
Figure 4: Cooling curves for a point at the center of a single bead deposit using several values of the spatial cutoff Rcut . The beam passes over the point in question at t0 . Example 3 decay profiles proportional to t−1 (dashed grey) and t− 2 (doted grey) are also indicated for reference.
d
tational expense with accuracy requirements. In the case of figure 4, Rcut =
te
1.0 mm ensures a reasonably accurate thermal history for elapsed times less than 20 ms, at which point the temperature rise induced by the source has dis-
Ac ce p
sipated to approximately 100 ◦C, well below both solid-liquid and solid-state transformations for Ti-6Al-4V. 3.3. Verification: Comparison to Rosenthal Solution The present model can be compared to analytic results for melt pool dimen-
sions in the case of a beam moving at uniform velocity. A series of single bead calculations with power ranging from 150 − 600 W, speed 500 − 2000 mm s−1 ,
and initial temperatures 25 − 1000 ◦C were executed with the DSM. The melt
pool widths were measured by finding points on the melting isotherm when the melt pool was in steady state, far from either the beginning or end of the track, and then compared to results predicted by equation 1. Figure 5(a) shows a comparison of the DSM results with the Rosenthal results and a linear least squares regression between the two. There is excellent agreement with the an-
18
Page 18 of 48
alytic result, with DSM widths generally within 1 µm (.1.0%) of the analytic results across the very wide range of P, v, and T0 considered. A linear regres-
ip t
sion between DSM and Rosenthal widths has slope of 0.996 and an intercept of only 0.35 µm. Similarly, figure 5(b) shows similar excellent agreement for
cr
the melt pool length calculation. Together, these results verify that when run with appropriate numeric parameters the present implementation of the DSM
us
reproduces expected analytic solutions.
For the single track configuration, one other main difference between the DSM and the Rosenthal solution is the finite size of the energy source as char-
an
acterized by the 4σ diameter. We define a P´eclet number for this problem 4σv 4σρcp v = . α κ
Pe =
(20)
M
Employing a 4σ beam diameter on the order of 100 µm and using thermophysical properties for Ti-6Al-4V we find Pe & 10 for v & 1 m s−1 . It is not until the velocity is decreased toward 0.1 m s−1 range that we would expect to see a
d
significant impact of finite spot size on melt pool dimensions.
te
3.4. Thermophysical properties
Another assumption in this model is the temperature independence of the
Ac ce p
parameters κ, ρ, and cp . As indicated in equations 1 and 2, w and l could be
affected by temperature dependence of the product ρcp , and κ respectively. Figure 6 shows the w and l of a melt pool measured at the center of the single bead
deposit described in Figure 1(a) as a function of the temperature Tp at which the fixed parameters are taken for Ti-6Al-4V. For this alloy, the resulting full variation in melt pool widths is on the order of 10% across the full range of property temperature. The length, however, is more sensitive, particularly when using property temperatures below 500 ◦C. This is, of course, alloy dependent, and for the duration of this work we will employ properties at Tp = 1203 ◦C. This value results in effective properties that match the total sensible heat required
19
Page 19 of 48
(a)
cr
700 600 500
us
400 300
1:1 DSM fit: y=0.996x + -0.355 DSM FDM Hsl = 0 FDM Hsl = 286kJ/kg
200 100 0
8000
200
300
400
500
Rosenthal Width [μm]
(b)
7000
600
700
800
d
6000 5000
te
4000 3000 2000
Ac ce p
Numeric Length [μm]
100
M
0
an
Numeric Width [μm]
ip t
800
1:1 DSM DSM fit: y=0.997x + -21.735 FDM Hsl = 0 FDM Hsl = 286kJ/kg
1000
0
0
1000
2000
3000
4000
5000
Rosenthal Length [μm]
6000
7000
8000
Figure 5: Single-bead steady-state melt pool width (a) and length (b) predicted by the DSM and the FDM (with and without Hsl ), vs. the Rosenthal solution[14]. Solutions computed 343 discrete combinations of absorbed power (P = 150 − 600 W), speed (v = 500 − 2000 mm s−1 ), and initial temperature (T 0 = 25 − 1000 ◦C) using thermal properties for Ti-6Al-4V. A linear least-squares regression indicates good agreement between the DSM and Rosenthal solution. The boxed region indicates the region applicable to the experimental results discussed in section 4.1. FDM results computed with Hsl = 0, and 286 kJ kg−1 . Realistic Hsl > 0 has a small impact on steady state width, but a significant impact on the observed melt pool length.
20
Page 20 of 48
3.5
290
3.0
285
2.5
280
2.0
275 270
1.0
Width Length 0
250
500
750
1000 1250 1500 1750 2000
Temperature of fixed parameters [C]
an
265
us
1.5
ip t
295
cr
4.0
Length [mm]
Width [μm]
300
M
Figure 6: Sensitivity of predicted single track melt pool width w and length l as a function of the temperature Tp chosen to evaluate the thermophysical parameters ρ, cp , κ, using parameters for Ti-6Al-4V from Ref. [36]. Vertical dotted lines at T = 1000 ◦C and 1650 ◦C are the approximate locations of the α ↔ β and β ↔ liquid phase transformations respectively.
d
to bring the material from T0 = 25 ◦C to the melting point, i.e.
te
ZTm ρ (Tp ) cp (Tp ) (Tm − T0 ) = ρ (T ) cp (T ) dT.
(21)
T0
Ac ce p
Fully temperature dependent thermophysical properties inherently make the governing equations non-linear, and thus not amenable to the solution approach utilized in the DSM. However, in most alloys this effect is of secondary importance and should not inhibit the use of the DSM for basic initial calculations. We already discussed the effect of Hsl via a simple argument based on the
St number of the alloy. Here, we consider a more complex approach to calculate the effect of the Hsl on aspects of the melt pool shape for the single track case, motivated by results in Section 3.2. It was seen there that for single tracks, diffusion along the direction of beam motion is suppressed, and the temperature initially decays as t−1 , effectlvely a 2D cooling problem. We formulate the diffusion problem in a cylindrical coordinate system (r0 , θ0 , z 0 ), with the cylinder axis z 0 parallel to the direction of beam travel. All derivatives in θ0 and z 0 are
21
Page 21 of 48
Width [µm] 211.1 209.7 211.6 204.8
Width Diff.[%] – -0.7 0.2 -3.0
Length [µm] 1238.9 1210.1 1242.1 1661.7
Length Diff.[%] – -2.3 0.3 34.1
ip t
Configuration Rosenthal DSM FDM, Hsl = 0 FDM, Hsl = 286 kJ kg−1
us
cr
Table 2: Comparison of length and width predicted by the Rosenthal model, the DSM, and the FDM formulated in cylindrical coordinates, in the latter case both with and without latent heat of fusion.
neglected, we set σy = σz , and solve the remaining 1-D diffusion equation with a simple finite difference model (FDM) formulation (see supplemental information
an
for further details). This method is strictly only applicable to sources moving in straight line at uniform velocity, and is included here merely to estimate the magnitude of the effect across a range of vayring P, v, T0 . The apparent heat
M
capacity methodology described in [44] and references therein is employed to account for non-zero Hsl . Table 2 summarizes the predicted steady state w and l of a single track melt pool using the standard conditions with the Rosenthal
d
solution, the DSM approach, and the FDM calculation with both Hsl = 0 as
te
well as a realistic Hsl = 286 kJ kg−1 .
As shown in table 2 and figure 5(a), with Hsl = 0 the FDM solution agrees
Ac ce p
well with the Rosenthal and DSM, the estimated width being generally within 1%, verifying the assumptions and implementation. Using a realistic Hsl = 286 kJ kg−1 in the FDM, the width is slightly decreased and the length and
duration molten grows by ≈ 34% compared to the Rosenthal and DSM lengths, as shown in figure 5(b). Note that in general the effect of Hsl will depend on the specific P, v, and T0 . In particular, neglecting Hsl is a poorer assumption as T0 becomes large because the effective St becomes smaller. For the remainder of the present work, Hsl is ignored. This likely has limited effect on the melt pool width for the alloy in question, but will lead to underestimating the length and melt-duration, on the order of 35%.
22
Page 22 of 48
3.5. Boundary Conditions The solutions in equations 14 and 19 are for a semi-infinite slab with uni-
ip t
form thermophysical parameters throughout the volume. The effective value of α in the powder bed surrounding a component is likely over-estimated in this
cr
configuration, leading to enhanced cooling particularly near the edges of a component. To assess the magnitude of this effect we consider the opposite limiting
us
case by including a zero-flux boundary condition, effectively enforcing zero thermal conductivity in the powder bed outside the part. There are various ways of incorporating boundary conditions into the DSM, but the most expedient for a
an
simple convex cross section is via fictitious “mirror” sources placed outside of the component.
Consider a square cross section 7 × 7 mm as shown in figure 1(c) with con-
M
ventional alternating or “snaking” scan vectors running parallel to x direction and spanning from edge to edge of the square. Figure 7 shows the thermal histories computed at two points on the top surface at positions x, y, z =
d
0.05, 3.50, 0.0 mm (just inside of the boundary), and x, y, z = 2.00, 3.50, 0.0 mm
te
(closer to the center of the feature). The thermal histories are computed with a zero-flux boundary condition applied 75 µm toward the powder-bed side of the
Ac ce p
component-bed boundary. This position is chosen to account for the fact that the melt pool typically extends beyond the position of the center of the beam; some over-melting and finite conductivity in the material immediately adjacent to the component will result. The specific position chosen here is approximately equal to half of a typical melt pool diameter observed for isolated single track cases. Solid lines in the figure are for the no-flux case, and dashed lines are for the original uniform diffusivity case with no additional boundary condition. The difference between the two conditions is shown as a dotted line on the right hand axis. For this analysis, the cutoff distance Rcut is set to a value much larger than the size of the feature such that all locations within the part could be influenced by the boundary. Table 3 summarizes multiple metrics of the difference between histories computed at a number of positions along the center-line of the square. Both the 23
Page 23 of 48
t∆Tmax − t0 [ms] 1.6 2.3 3.2 19.6 53.3 315.8 606.2 722.7
R
∆T dt[◦C s] 248 245 242 236 227 215 202 196
ip t
T∆Tmax [◦C] 1594 1392 1198 662 440 205 151 140
cr
∆Tmax [◦C] 461 341 268 164 105 65 54 51
us
x[mm] 0.00 0.05 0.10 0.25 0.50 1.00 2.00 3.50
an
Table 3: Metrics computed for thermal history differences at several positions x from the boundary for the square geometry described in figure 7, all at y = 3.50 mm). ∆Tmax gives the largest instantaneous temperature difference between the case with and without the no-flux boundary condition. T∆Tmax is the temperature of the no-flux condition solution when the maximum error occurred, t∆Tmax − t0 is elapsed time between the source directly interacting with the point in question and the time when the maximum discrepancy is observed, and R ∆T dt is the integral of the temperature discrepancy over the full history.
M
table and figure 7 indicate that that for a point near the edge, the effect of the boundary condition is generally stronger and is felt at a time closer to the primary energy input event. The maximum difference (∆Tmax = L∞ norm) be-
d
tween the temperature histories is 461 ◦C on the edge of the square (x = 0.0 mm)
te
and decreases toward the center to 51 ◦C. The temperatures T∆Tmax and times t∆Tmax − t0 of the no-flux boundary condition solution at which these maximum
Ac ce p
differences are observed are also tabulated as a function of distance from the edge. It is clear that T∆Tmax decreases as the distance of the position in question from the edge increases, and for distance & 500 µm from the edge, T∆Tmax is
low enough that the effect is likely negligible for this alloy. For the duration of this work we will neglect the effect of the boundary conditions except where stated otherwise, though acknowledge that their influence affects the buildup of heat over long times and their significance is certainly a function of component geometry.
3.6. Time cutoff In Section 3.2 we considered the effect of a varying spatial cutoff Rcut with no time cutoff, tcut = ∞. Here we consider the opposite limiting case with varying tcut and Rcut = ∞. Figure 8 shows the effect of tcut on the thermal
24
Page 24 of 48
ip t cr
us 2500 2000
300
250
T, Temperature [C]
an
T [C]
103
350 T, No flux, x= 0.05[mm] T, No BC, x= 0.05[mm] ΔT, No flux, x = 0.05[mm] T, No flux, x= 2.00[mm] T, No BC, x= 2.00[mm] ΔT, No flux, x = 2.00[mm]
1500
M
1000
500
150 0.150
0.155
0.160
0.165
Time [s]
0.170
0.175
0.180 100
d
102
200
te
50
Ac ce p
0
10−1
100
101
Time [s]
Figure 7: Thermal histories computed at a position near the edge (x, y = 0.05, 3.50 mm, blue) and closer to the center (x, y = 2.00, 3.50 mm, red) of a 7 mm square with a no-flux boundary condition (solid) applied at an offset of 75 µm outside the square, or no boundary condition (dashed) (e.g. spatially constant thermophysical properties in the surrounding bed). Dotted lines indicated the temperature discrepancy, ∆T , between the two solutions on the right hand axis. Inset: zoom of thermal histories in a time window close to the primary energy input event for the positions in question. Note that no distance cutoff is used for this calculation.
25
Page 25 of 48
ΔT, Temperature Error [C]
3000
history for the point at the center of the square feature in figure 1(c), both for the full time history in figure 8(a) and in the region near the primary energy
ip t
input events in figure 8(b). tcut = ∞ is a bounding case with no time cutoff, and the histories increasingly diverge from this case as tcut decreases.
cr
Comparing figure 8 to 4, we can see that finite tcut induces a more rapid error accumulation than Rcut . In the case of finite Rcut , the divergence of the thermal
us
history from the limiting condition at a point ~rj is driven by artificial removal of energy at locations distant from ~rj . This in turn enhances the effective rate of cooling at ~rj , which will, in the limit, lead to cooling that follows t−3/2 behavior.
an
However, in the case of tcut the divergence is a result of artificial energy removal in the immediate vicinity of ~rj , and thus the effect is felt essentially immediately. The rate of error accumulation due to tcut is therefore controlled by the geometry
M
and sequencing of the particular scan strategy employed, not diffusion. Additionally, while the effect of removing an individual source might be small based on the analysis in Section 2.5, if a particular scan strategy results in many
d
points that are spatially and temporally close to each other, their individual
te
weak impacts can still add into large temperature effect. In practice, finite values of both tcut and Rcut are employed, and thus a source could be excluded
Ac ce p
by either threshold. For the reasons noted above, specific choices of tcut and Rcut must be considered jointly, and their their impacts assessed for different scan strategies. In the present case with a single, continuous, and smoothly moving source, the choice tcut = 0.1 s is effective at tracking the history to
temperatures below 500 ◦C, well below the melting point and important solidstate transformations for Ti-6Al-4V. 4. Experimental Measurements Experiments were executed to generate calibration and validation data for the model described above. Parallel single beads were deposited under a range of nine unique P, v combinations, with 2 additional replicates at a centerpoint condition, for a total of 11 beads. An EOS M280 laser PBF machine melted individual vectors after a Ti-6Al-4V powder layer of nominal thickness 30 µm 26
Page 26 of 48
10−1
t [s]
te
1400 1200 1000
Ac ce p
T [C]
100
d
(b)
M
an
102
1600
ip t
us
tcut = ∞ tcut = 1.0 tcut = 0.3 tcut = 0.1 tcut = 0.03
cr
(a)
T [C]
103
tcut = ∞ tcut = 1.0 tcut = 0.3 tcut = 0.1 tcut = 0.03
800 600 400
0.14
0.15
0.16
0.17 t [s]
0.18
0.19
0.20
Figure 8: Effect of varying tcut on the thermal history of the central point in figure 1(c) over the full duration of the scan (a) and a zoom near the most direct energy input event (b). Curves for tcut = ∞, 3.0 s and 1.0 s are effectively equivalent within the window shown in (b).
27
Page 27 of 48
was spread onto a wrought Ti-6Al-4V substrate, initially at room temperature. The individual beads were 40 mm in length and separated by 8 mm. The 4σ
ip t
spot diameter was approximately 75 µm as characterized by beam profiling, and the beam trajectory and speed were verified by an in situ monitoring system[45].
cr
The typical time elapsed between the start of adjacent beads is 0.14 s, and the
total time to deposit all eleven beads was 1.54 s. For comparison, the diffusion
us
time scale l2 / (4α) assuming properties at room temperature for this alloy across a distance l = 8 mm is 2.35 s, and thus the beads are effectively thermally independent on the timescale of the deposition.
an
The resulting single beads were characterized via optical and scanning electron microscopy. Top-down Backscatter Electron (BSE) image montages of single beads were collected in a Tescan Lyra dual-beam FIB. Bead width transients
M
were observed to extend typically no more than ≈ 1mm from the beginning of the bead, and thus the mean and standard deviation of bead width (w, ¯ σw ) were assessed across a steady-state region of interest between 2 and 8 mm from the
d
beginning of the beads. The width was recorded at approximately 30 discrete
te
locations each separated by 200 µm.
Additionally, approximately 15 cross-sections were collected within the same
Ac ce p
steady-state region for each bead via serial sectioning [46]. Cross-sections were etched with Kroll’s reagent (92% H2 O, 6% nitric acid, 2% hydroflouric acid) and optically imaged to reveal the melt pool and surrounding heat affected zone. The mean and standard deviation of the width (w, ¯ σw ), depth of the ¯ σd ), and height above minimum melted location below the substrate surface (d, ¯ σh ) were characterized from these images by three independent the substrate (h, observers, and the results were averaged together. Several examples of crosssection images are shown in the panels of figure 10. The typical variation in the mean width as determined via the top-down compared to cross sections was approximately 1%, and for no bead was the variation > 3.5%. Figure 9 shows a comparison of the top down and cross section measurements, along with bars indicating the measurement standard deviations.
28
Page 28 of 48
180
y
= 0.966x + 3.231
ip t
160
cr
150 140 130 120
130
140
150
160
Cross Sectional Width [μm]
170
180
an
120
us
Top Down Width [μm]
170
M
Figure 9: Bead width as measured via cross sections and top down views. Bars indicate 1 standard deviation on either side of the mean value.
4.1. Single Track Width and Depth
d
Experimentally measured width and depth are shown in Table 4. As ex-
te
pected, both width and depth increase as either power is increased or speed is decreased. Additionally, the depth to width ratio ranges from 0.678 to 0.921 (a
Ac ce p
circular cross section is characterized by d/w = 0.5). In general, the depth to width ratio increases with increasing power, or decreasing speed. The value of the source depth to radius utilized in the model, σz /σx , must be calibrated to
match the observed depth to width ratio. The choice σz /σx = 1 results in melt
pools with circular cross sections and thus σz /σx must be increased to match
the observed melt pool cross section behavior. Furthermore, experimental data indicates considerable variation in w at a
fixed P/v. In particular, 3 replicate tracks at the center-point P, v conditions were produced, and observed mean widths are shown in Table 4, as condition IDs 1, 2, and 9, annotated with asterisks. Tracks 1 and 2 demonstrated similar widths, whereas track 9 exhibited a width ≈ 10% lower. For the purposes of calibrating the model and computing error between empirical observations and model outputs, we utilize the values from track ID 2 for this particular P, v 29
Page 29 of 48
ip t cr us an M d te Ac ce p
Figure 10: Examples of single bead cross sections utilized to determine width and depth. (a) condition 11 (218 W, 975 mm s−1 ), (b) condition 2 (290 W, 1300 mm s−1 ), and (c) condition 10 (363 W, 1625 mm s−1 ). Approximately 30 top down measurements and 15 cross section measurements were collected for each P, v condition.
30
Page 30 of 48
Depth[µm] d¯ σd 100.4 9.2 112.1 4.3 95.3 3.1 141.5 1.8 90.0 4.2 133.9 5.8 86.2 3.7 150.2 3.9 116.9 2.8 121.2 3.7 109.3 3.1
Height[µm] ¯ h σh 25.9 5.7 21.7 3.4 29.4 24.0 20.7 3.2 19.0 6.7 22.1 1.9 26.4 9.8 16.0 4.0 13.6 4.6 22.1 11.1 11.5 2.0
ip t
Width: Top[µm] w ¯ σw 144.1 6.5 143.3 7.1 125.4 8.4 155.4 6.3 125.4 5.3 175.2 5.4 124.4 6.4 172.0 6.4 129.6 4.8 130.6 5.5 152.0 3.8
cr
Width: Sec.[µm] w ¯ σw 147.9 7.0 144.9 5.2 127.2 5.6 153.9 4.2 126.0 5.4 177.3 6.7 122.2 4.8 174.2 8.0 134.0 4.8 131.6 3.5 156.2 4.7
v[mm s−1 ] 1300 1300 1860 1300 1300 908 1733 1040 1300 1625 975
us
P [W] 290 290 341 341 239 239 290 290 290 363 218
an
ID 1* 2* 3 4 5 6 7 8 9* 10 11
M
Table 4: Mean and standard deviation across all measurements of the width (w, ¯ σw ), depth ¯ σd , and height h, ¯ σh of single tracks 1 through 11, in micrometers. Width is characterd, ized by via both cross section (Sec.), and top-down images (Top), whereas height and depth are determined from cross section images only. * denotes replicates at the center-point of the processing conditions.
combination, as it is the most regular of the three deposits. Further work is
d
needed to understand the source of bead-to-bead variation in cases with nominally identical P and v values, but we note that these differences are still small
te
compared to the overall range of widths observed across the full range of P and v in the experiments.
Ac ce p
4.2. Model Calibration
As described in Sec. 2.1, two adjustable parameters are available to bring
the predicted melt pool geometry into agreement with observations. Specifically, the overall efficiency factor η and the volumetric source’s depth to radius ratio
σz /σx must be adjusted to reproduce single track width and depth observations as described in table 4.
Calculations of the melt pool width and depth were performed with the
DSM using the same geometry and processing conditions as the experimental single beads described above. The melt-pool dimensions are determined from extents of the isotherm at Tm as in section 3.3. A range of η and σz /σx values were employed, and width predictions resulting from optimal choices of these calibration parameters are presented in figure 11 along with experimentally 31
Page 31 of 48
Experimental
180
cr
Width
[μm]
200
σz/σx only η only σz/σx and η
ip t
220
160
120
0.18
0.20
0.22 0.24 P/v [J/mm]
0.26
0.28
0.30
an
0.16
us
140
M
Figure 11: Single bead measurements and calculations with varying degrees of calibration. Linear fits over a narrow range of P/v are included as a guide to the eye.
observed widths as a function of linear energy density P/v. Note that P is the
d
nominal power setting before any efficiency factor is applied.
te
Figure 11 shows the model performance under three different calibration cases. In the first case σz /σx is adjusted at constant η = 1 until the depth/width
Ac ce p
ratio best matches observations at the center-point in the P, v space. In the second case, η is tuned to match the observed width with fixed σz /σx = 1.0. In the final case, σz /σx is first adjusted until the depth to width ratio matches observations at the center-point in the P, v space and then η is adjusted until the width matches across the P/v range. The best agreement between modeled and observed width occurs when both
η and σz /σx are tuned, as indicated in figure 12. It should be noted that the optimal choice of σz /σx varies as a function of P and v, but here we only optimize for the centerpoint of the P, v values. Note that equation 1 indicates that w scales as
p P/v. However, the range of
P/v in figure 11 is relatively narrow and does not exhibit extensive curvature and thus linear fits are used to characterize the relationship. Care must be exercised if extrapolating beyond the range of observed data, though calculations with 32
Page 32 of 48
190
ip t
180
160
140
cr
150
us
/
σz σx only η only σz σx and η
130
/ 1:1 1:1±5,10%
120 110 110
120
130
140
150
160
Experimental Width [μm]
170
180
an
Model Width [μm]
170
M
Figure 12: Experimental vs. modeled widths for single tracks across all power, speed conditions for several values the efficiency factor η, and source depth to radius ratio σz /σx .
d
the DSM over a much wider range of P and v values have been executed and p confirm the expected P/v scaling.
te
As expected, the model over-predicts the bead width significantly with η = 1.0. Decreasing the efficiency while keeping σz = σx brings the predictions in
Ac ce p
reasonably good agreement with the observed widths, though the slope of the width with respect to P/v is underestimated. Finally, tuning both η and the σz /σx ratio to match the width at the center-point and brings the slope of the w vs. P/v into better agreement with observations. Again, we point out that previous experimental work has indicated that
the absorption efficiency can be a strong function of the processing conditions, particularly as power is increased near the onset of keyhole instability [34]. Thus, we do not expect a single combination of η, σz /σx to work over a wide range of processing conditions. Instead, the model must be calibrated over the full intended use range of P and v. For the remainder of this work, we will
employ P and v conditions representative of the centerpoint in the range shown in figures 11 and 12, and thus will utilize the η and σz /σx tuned to the center point of this range. 33
Page 33 of 48
5. Model Validation
ip t
After calibrating the DSM against single track experimental results, it is tested against empirical observations of melt pool geometry on the top surface
of an object produced by a multi-vector scan pattern. The test geometry is a
cr
rectangular prism with dimensions of 40 × 10 mm in the build plane, and 20 mm
along the build direction, shown in figure 13. The scan pattern used to process
us
the topmost layer of the item is shown in figure 13(a), colored by the position of the laser as a function of time. The top surface is divided into sub-domains known as ‘stripes’ of width 7 mm. The stripe boundaries are inclined to the long
an
axis of the component by 68.26◦ , and are clearly visible in the BSE montage of the object as a series of dark lines crossing the item in figure 13(b). The
M
individual scan vectors within the boxed region of interest are orthogonal to the stripe boundaries, and are shown in more detail in figure 14(a). When the beam reaches the end of each vector it is shut off before turning and accelerating to
d
begin scanning the next vector, at which point it is turned back on. The mean empirically observed turning time is employed for all turning events in the path,
te
and the beam is assumed to move at the nominal processing speed as it traverses from the final vector in a stripe to the first vector in the subsequent stripe.
Ac ce p
A calculation with this scan pattern is executed with the parameters listed
in table 1 with the exception of the following modifications; the spatial cutoff is set to Rcut = 4 mm to ensure the boundary condition can affect internal regions consistent with sections 3.2 and 3.5, tcut = 0.1 s as described in section 3.6, and T0 = 100 ◦C which is representative of a multilayer deposit. A no-flux boundary condition is employed at a distance 150 µm outside all free edges of the item. The specific region of interest shown in figure 14 is selected because the scan-
ning vectors interact with both a free surface (the lower edge), as well as the internal stripe boundary, annotated as a solid black line labeled “0”. The vectors are sequentially scanned beginning in the upper right of the region and moving toward the lower left. While the processing velocity is nominally uniform, because the vectors become successively shorter as the corner is approached the
34
Page 34 of 48
ip t cr us an
M
Figure 13: (a) Global scan path covering a 40 × 10 mm rectangular prism with scan vectors inclined 21.74◦ from the horizontal. Color indicates the time elapsed since the beginning of the scan on the layer. (b) BSE image montage covering the full object. Boxes in (a) and (b) indicate the region of interest shown in Figure 14, where individual scan vectors are shown.
d
time between reheating events becomes shorter, and heat begins to build up.
te
Well away from either the internal stripe boundary or the free edge of the item, measurements from figure 14(a) indicate that the un-remelted width of the pool perpendicular to the scan vectors achieves a steady state value
Ac ce p
138.1 ± 10.0 µm, approximately equal to the expected nominal spacing between
the parallel scanning vectors of 140 µm, which is also the value that would be predicted by the Rosenthal model. However, near both the internal stripe boundary and the free surface, the melt pool width is larger for vectors which are moving away from the boundaries, sometimes by a significant amount. This locally increased width in turn partially covers the adjacent, previously processed vector making the apparent width of vectors approaching the boundary appear narrower than nominal, creating a pattern of alternating wide and narrow tracks. Figure 14(b) shows model output, with predicted thermal histories postprocessed to identify the un-remelted width of melt pools on the top surface at
35
Page 35 of 48
ip t cr us an M d te Ac ce p
Figure 14: Region of interest identified in figure 13. (a) BSE image showing spatial variation in the un-remelted apparent width of the scan vectors, along with scan vectors colored by time processed, and measurement location lines at 0.3, 1.5, and 3.0 mm [red, green, blue respectively] from the stripe boundary [black]. (b) DSM prediction of un-remelted top surface with the same vectors and measurement lines. Larger versions are available in the supplemental information.
36
Page 36 of 48
Experiment DSM
200 150 50
−100
Line 1: 0.3mm
−150 100
us
50 0 −50 −100
Line 2: 1.5mm
an
−150 50 0 −50 −100
Line 3: 3.0mm
−150 5
10
Scan Vector
15
20
te
d
0
M
Deviation from Ideal Unremelted Width [μm]
0 −50
cr
ip t
100
Ac ce p
Figure 15: Measurements of the un-remelted top surface in the region shown in figure 14 as predicted by the DSM (open symbols) and measured from the experiment (closed circles), expressed in terms of deviation from the hatch spacing, 140 µm, which is the prediction of the Rosenthal solution. The grey band indicates typical variation in observed width well away from a stripe boundary or free surface of ±10.0 µm.
z = 0, and figure 15 shows the un-remelted width measurements from both the experimental image as well as the model output along measurement lines labeled 1,2,3. The scan vectors have been discretized in the same fashion as employed in all earlier examples. The width is expressed in terms of deviation from the nominal 140 µm hatch spacing for each scan vector crossed by the measurement
line, beginning with vector 1 nearest to the free surface, and counting into the part. For Line 1 [red] at 300 µm from the stripe boundary, both the observed and modeled widths differ significantly from the nominal 140 µm hatch spacing, systematically alternating between narrower and wider. As described above, 37
Page 37 of 48
vectors moving away from the stripe boundary experience a pre-heat effect due to the previously scanned vector. They partially, or in some cases completely,
ip t
cover the previously deposited adjacent track, leading to the alternating wide-
narrow-wide behavior. Measurements in figure 15 show that the DSM also
cr
captures this trend, though the amplitude of the oscillation is under-predicted. Additionally, the vectors closest to the free surface on Line 1 show that as
us
the corner is approached, the overall scan vector length decreases, decreasing the time between re-heating events, increasing preheat, and therefore width. Qualitatively, it can also be seen that both the experiment and DSM predict
an
temporal overlap of melt pools from adjacent vectors, for the vectors closest to the intersection of the stripe-boundary and free surface. Further from the
as the beam is turning around.
M
corner, the melt pool induced by a particular track is able to completely solidify
Similar trends are observed for measurements along Line 2 [green] at 1.0 mm from the stripe boundary, though the amplitude of the alternating width is
d
weaker, as the pre-heat variation is generally smaller, again, with the DSM
te
under predicting the variation. Finally, for Line 3 [blue] at 3.0 mm from the stripe boundary, the observed and predicted width are both generally close to
Ac ce p
the expected value of 140 µm, with deviations apparent only for scan vectors close to the free-surface. Variations in the width observed along Line 3 in the DSM result are generally smaller than the experimentally observed variation in track width of ±10 µm in regions far from the stripe boundary or free surface. While the DSM is able to capture a number of trends in the spatial variation
of melt pool shape and duration molten, there are still discrepancies with the empirical results. Fluid flow, particularly near the beginning of a scan vector, is likely important in determining the extent to which the pool flows over the top of the adjacent track, with the expectation that flow driven by recoil pressure would drive liquid to recover adjacent tracks to a higher degree than predicted by the DSM. Latent heat of fusion was observed to have a slight impact on single track widths, but a more significant effect on pool lengh, and therefore duration molten. Longer lived molten pools would have more of an opportu38
Page 38 of 48
nity to overlap with those from adjacent tracks, again leading to a larger molten pool than predicted by the DSM. Additionally, the actual boundary condition is
ip t
somewhere between zero flux and uniform conductivity. However, similar calculations without the no-flux condition do not show a noted change in unremelted
cr
width for the end of vector features along the stripe-boundary, and thus is likely not the primary cause of discrepancy. Finally, we employ a uniform value of
us
T0 = 100 ◦C. We expect in practice that this value will have some spatial variation, however the inter-layer times are significant and the geometry of this build
an
simple enough that it is unlikely to drive the observed discrepancies. 6. Performance and Scaling
The DSM has been implemented as a filter within the DREAM.3D software
M
environment[42]. Specifically, the model has been implemented in C++ with thread parallelization performed using Intel’s Threaded Building Blocks (TBB). In an effort to quantify the practical scalability of the model, we performed
d
both a strong and weak scaling study. The study was carried out on an HP
te
Z840 with dual Intel Xeon E5-2687W V4 CPUs with base frequency 3.0 GHz, resulting in 48 logical processors. The strong scaling study involved evaluating
Ac ce p
3.75M spatial locations utilizing nLP =1, 6, 12, 24 and 48 logical processors. The weak scaling study utilized the same nLP values, but adjusted the number
of spatial evaluation locations Nj between 0.3125M, 1.875M, 3.75M, 7.5M and 15M correspondingly. The results of the two studies are listed in table 5, both in actual wall time and in scaling efficiency. For the weak scaling study, perfect load balancing was not explicitly enforced, rather TBB’s default load balancing scheme was used. The specific problem geometry and model parameters used are the same as described in section 5. The time resolution of the calculated thermal histories at each spatial evaluation location was selected to be 100 µs in all cases. As described in Section 2.5, the source discretization and output time resolution are held constant, and thus each spatial output location has, on average, the same number of evaluations of equation 9.
39
Page 39 of 48
Both strong and weak scaling efficiencies indicate that the DSM scales well. In both cases, scaling efficiency of ≈40% is achieved for the extreme cases. Fur-
ip t
thermore, while untested here, efficient scaling of parallelization in distributed
memory environments is also expected due to the uncoupled nature of the model,
cr
and is a means to further increase feasible values of Nj .
The case with Nj = 15 × 106 and nLP = 48 evaluated approximately 6300
us
output points per second of wall time on the hardware employed. Accounting for a 3 dimensional distribution of source points at the same density employed for the calculation in figure 14, using Rcut = 1 mm, the DSM could screen
an
temperature histories at points separated by 1 mm filling a cube with edge length 25 mm in approximately 60 s of wall time. Alternately, filling a cube 250 mm on edge (approximately a full build volume) with the same 1 mm grid of points
M
would require 15 h. For comparison, we estimate that an actual build of this height would require approximately 140 h of printing time using a generous estimate of 60 s per layer. Additionally, we note that generally significantly
d
less than 50% of build volumes are typically occupied with parts, and thus the
above.
te
spatial resolution could be increased while maintaining the same wall time cost
Ac ce p
There are opportunities to sample the scan vectors and place source points more efficiently than the uniformly spaced approach taken here with little or no loss of accuracy. Future work will focus on characterizing and quantifying such schemes, but initial estimates indicate a further reduction in the above mentioned 15 h calculation to approximately 2 h. Finally, while the above examples place the evaluation points on a cubic grid for simple comparison, we note again that there is complete freedom in placement and number of these points within the volume of interest, with no loss of accuracy for any particular thermal history. 7. Conclusion A discrete source model for powder bed fusion thermal history has been developed and verified against analytic calculations for single track geometries. 40
Page 40 of 48
6 twall E 1121 94 2167 95 -
12 twall 1179 -
E 89, 87 -
24 twall 971 1712 -
E 53 62 -
48 twall 655 2365
E 39 45
ip t
1 twall E 1053 100 12378 100 -
cr
nLP → N j × 10−6 0.3125 1.875 3.75 7.5 15
us
Table 5: Wall Time (twall ) and scaling efficiency (E) results for a strong and weak scaling study using several numbers of logical processors nLP . Wall Times are in seconds, efficiencies are in percents, and numbers of spatial evaluation points Nj are in millions. Note the two efficiencies listed for the nLP = 12 case are for the weak and strong scaling, respectively.
an
Convergence studies have been presented, and guidelines for choosing parameters to ensure required levels of accuracy have been developed, and an understanding of the range of approximation errors have been quantified.
M
A model calibration methodology using the geometry of experimental singletrack deposits has been presented, and it was shown that the DSM can be calibrated to match both the melt pool width and depth for Ti-6Al-4V.
d
The DSM was utilized to predict several spatial and temporal transient ef-
te
fects due to thermal interactions induced by multi-track deposits. Scan pattern and component geometry combine to produce spatial variation in pre-heat, and
Ac ce p
thus track geometry. Several reasons for quantitative discrepancies between experiment and DSM predictions have been presented. The DSM is an effective means of quickly generating an estimate of the local
thermal history induced by complex scan strategies when combined with arbitrary component geometry. We expect this tool will be useful for screening and identification of situations leading to anomalous thermal history, and ultimately to direct usage of higher fidelity models and experimental capabilities. Simple temporal discretization schemes for the energy source have been utilized here for clarity, but more optimized sampling for generation of source points in part scale calculations is expected to further enhance performance and is currently being pursued. Additionally, the thermal histories resulting from the present model are still cumbersome to work with, and reduced dimensional representations of
41
Page 41 of 48
thermal history will be a thrust of future investigation.
ip t
Acknowledgments
We wish to acknowledge Jonathan Miller at the Air Force Research Labo-
cr
ratory for many helpful conversations. Thomas Spears of GE Aviation and Joy
Gockel of Wright State University were instrumental in executing single track
us
deposition experiments, and Sonya Boone of UES collected images of multi-track deposits.
an
References
[1] T. J. Horn, O. L. Harrysson, Overview of current additive manufacturing technologies and selected applications, Science progress 95 (3) (2012) 255–
M
282.
[2] W. E. Frazier, Metal Additive Manufacturing: A Review, Journal of
d
Materials Engineering and Performance 23 (6) (2014) 1917–1928. doi:
te
{10.1007/s11665-014-0958-z}.
[3] Y. Huang, M. C. Leu, J. Mazumder, A. Donmez, Additive Manufacturing:
Ac ce p
Current State, Future Potential, Gaps and Needs, and Recommendations, Journal of Manufacturing Science and Engineering - Transactions of the ASME 137 (1). doi:{10.1115/1.4028725}.
[4] W. E. King, A. T. Anderson, R. M. Ferencz, N. E. Hodge, C. Kamath, S. A. Khairallah, A. M. Rubenchik, Laser powder bed fusion additive manufacturing of metals; physics, computational, and materials challenges, Applied Physics Reviews 2 (4). doi:{10.1063/1.4937809}.
[5] R. Jiang, R. Kleer, F. T. Piller, Predicting the future of additive manufacturing: A delphi study on economic and societal implications of 3d printing for 2030, Technological Forecasting and Social Change 117 (Supplement C) (2017) 84 – 97. doi:10.1016/j.techfore.2017.01.006.
42
Page 42 of 48
[6] W. King, A. T. Anderson, R. M. Ferencz, N. E. Hodge, C. Kamath, S. A. Khairallah, Overview of modelling and simulation of metal powder bed
ip t
fusion process at Lawrence Livermore National Laboratory, Materials Science and Technology 31 (8) (2015) 957–968. doi:10.1179/1743284714Y.
cr
0000000728.
[7] M. J. Matthews, G. Guss, S. A. Khairallah, A. M. Rubenchik, P. J. Depond,
us
W. E. King, Denudation of metal powder layers in laser powder bed fusion processes, Acta Materialia 114 (2016) 33–42. doi:{10.1016/j.actamat.
an
2016.05.017}.
[8] S. Ly, A. M. Rubenchik, S. A. Khairallah, G. Guss, M. J. Matthews, Metal vapor micro-jet controls material redistribution in laser powder bed fusion
s41598-017-04237-z.
M
additive manufacturing, Scientific Reports 7 (2017) 4085. doi:10.1038/
[9] S. A. Khairallah, A. Anderson, Mesoscopic simulation model of selec-
d
tive laser melting of stainless steel powder, Journal of Materials Process-
2014.06.001.
te
ing Technology 214 (11) (2014) 2627–2636. doi:10.1016/j.jmatprotec.
Ac ce p
[10] S. A. Khairallah, A. Anderson, A. M. Rubenchik, J. Florando, S. Wu, H. Lowdermilk, Simulation of the main physical processes in remote laser penetration with large laser spot size, AIP Advances 5 (4) (2015) 047120. doi:10.1063/1.4918284.
[11] S. A. Khairallah, A. T. Anderson, A. Rubenchik, W. E. King, Laser powderbed fusion additive manufacturing: Physics of complex melt flow and formation mechanisms of pores, spatter, and denudation zones, Acta Materialia 108 (2016) 36–45. doi:{10.1016/j.actamat.2016.02.014}.
[12] B. A. Cowles, D. G. Backman, R. E. Dutton, Update to recommended best practice for verification and validation of ICME methods and models for aerospace applications, Integrating Materials and Manufacturing Innovation 4 (1). doi:10.1186/s40192-014-0030-8. 43
Page 43 of 48
[13] D. Rosenthal, The theory of moving sources of heat and its application to
ip t
metal treatments, ASME, 1946. [14] M. Tang, P. C. Pistorius, J. L. Beuth, Prediction of lack-of-fusion porosity
39 – 48. doi:10.1016/j.addma.2016.12.001.
cr
for powder bed fusion, Additive Manufacturing 14 (Supplement C) (2017)
[15] J. Raplee, A. Plotkowski, M. M. Kirka, R. Dinwiddie, A. Okello, R. R.
us
Dehoff, S. S. Babu, Thermographic Microstructure Monitoring in Electron Beam Additive Manufacturing, Scientific Reports 7 (2017) 43554. doi:
an
10.1038/srep43554.
[16] S. P. Narra, R. Cunningham, J. Beuth, A. D. Rollett, Location specific solidification microstructure control in electron beam melting of ti-6al-4v, Ad-
M
ditive Manufacturing 19 (2018) 160–166. doi:{10.1016/j.addma.2017. 10.003}.
d
[17] E. R. Denlinger, J. C. Heigel, P. Michaleris, T. A. Palmer, Effect of interlayer dwell time on distortion and residual stress in additive manufacturing
te
of titanium and nickel alloys, Journal of Materials Processing Technology
Ac ce p
215 (2015) 123–131. doi:{10.1016/j.jmatprotec.2014.07.030}. [18] T. Eagar, N. Tsai, Temperature-Fields Produced by Traveling Distributed Heat-Sources, Welding Journal 62 (12) (1983) S346–S355.
[19] J. Goldak, A. Chakravarti, M. Bibby, A new finite element model for welding heat sources, Metallurgical Transactions B 15 (2) (1984) 299–305. doi:10.1007/BF02667333.
[20] N. Nguyen, A. Ohta, K. Matsuoka, N. Suzuki, Y. Maeda, Analytical solutions for transient temperature of semi-infinite body subjected to 3-D moving heat sources, Welding Journal 78 (8) (1999) 265S–274S. [21] N. Nguyen, Y. Mai, S. Simpson, A. Ohta, Analytical approximate solution for double ellipsoidal heat source in finite thick plate, Welding Journal 83 (3) (2004) 82S–93S. 44
Page 44 of 48
[22] M. Van Elsen, M. Baelmans, P. Mercelis, J. P. Kruth, Solutions for modelling moving heat sources in a semi-infinite medium and applications to
ip t
laser material processing, International Journal of Heat and Mass Trans-
fer 50 (23-24) (2007) 4872–4882. doi:{10.1016/j.ijheatmasstransfer.
cr
2007.02.044}.
[23] M. Akbari, D. Sinton, M. Bahrami, Moving heat sources in a half space: Ef-
us
fect of source geometry, Proceedings of the ASME 2009 Heat Transfer Summer Conference HT2009 (2009) 685–694doi:{10.1115/HT2009-88562}.
an
[24] R. Komanduri, Z. B. Hou, Thermal analysis of the arc welding process: Part i. general solutions, Metallurgical and Materials Transactions B 31 (6) (2000) 1353–1370. doi:10.1007/s11663-000-0022-2.
M
[25] R. Komanduri, Z. Hou, Thermal analysis of the laser surface transformation hardening process, International Journal of Heat and Mass Transfer 44 (15)
d
(2001) 2845–2862. doi:10.1016/S0017-9310(00)00316-1. [26] R. Komanduri, Z. B. Hou, Thermal analysis of the arc welding process:
te
Part ii. effect of variation of thermophysical properties with temperature, Metallurgical and Materials Transactions B 32 (3) (2001) 483–499. doi:
Ac ce p
10.1007/s11663-001-0034-6.
[27] J. Winczek, Analytical solution to transient temperature field in a halfinfinite body caused by moving volumetric heat source, International Journal of Heat and Mass Transfer 53 (25-26) (2010) 5774–5781. doi: {10.1016/j.ijheatmasstransfer.2010.07.065}.
[28] S. Mohanty, J. Hattel, Cellular scanning strategy for selective laser melting Capturing thermal trends with a low-fidelity, pseudo-analytical model, Mathematical Problems in Engineering 2014 (2014) 715058. doi:10.1155/ 2014/715058. [29] S. Mohanty, J. H. Hattel, Cellular scanning strategy for selective laser melting: Generating reliable, optimized scanning paths and processing param45
Page 45 of 48
eters, Proc. SPIE Volume 9353, Laser 3D Manufacturing II 9353 (2015)
ip t
93530U. doi:10.1117/12.2079957. [30] S. Mohanty, J. H. Hattel, Reducing residual stresses and deformations in
selective laser melting through multi-level multi-scale optimization of cel-
III 9738 (2016) 93780Z. doi:10.1117/12.2212490.
cr
lular scanning strategy, Proc. SPIE Volume 9738, Laser 3D Manufacturing
us
[31] A. Plotkowski, M. Kirka, S. Babu, Verification and validation of a rapid heat transfer calculation methodology for transient melt pool solidification
an
conditions in powder bed metal additive manufacturing, Additive Manufacturing 18 (Supplement C) (2017) 256 – 268. doi:10.1016/j.addma. 2017.10.017.
M
[32] C. D. Boley, S. A. Khairallah, A. M. Rubenchik, Calculation of laser absorption by metal powders in additive manufacturing, Applied Optics 54 (9)
d
(2015) 2477–2482. doi:10.1364/AO.54.002477. [33] A. Rubenchik, S. Wu, S. Mitchell, I. Golosker, M. Leblanc, N. Peterson,
te
Direct measurements of temperature-dependent laser absorptivity of metal powders, Applied Optics 54 (24) (2015) 7230–7233. doi:10.1364/AO.54.
Ac ce p
007230.
[34] J. Trapp, A. M. Rubenchik, G. Guss, M. J. Matthews, In situ absorptivity measurements of metallic powders during laser powder-bed fusion additive manufacturing, Applied Materials Today 9 (Supplement C) (2017) 341 – 349. doi:10.1016/j.apmt.2017.08.006.
[35] S. Lu, H. Fujii, K. Nogi, Sensitivity of Marangoni convection and weld shape variations to welding parameters in O-2-Ar shielded GTA welding, Scripta Materialia 51 (3) (2004) 271–277. doi:10.1016/j.scriptamat. 2004.03.004.
[36] K. C. Mills, Recommended Properties of Commercial Alloys, Woodhead Publishing Ltd. & ASM International, 2002. 46
Page 46 of 48
[37] R. R. Dehoff, M. M. Kirka, W. J. Sames, H. Bilheux, A. S. Tremsin, L. E. Lowe, S. S. Babu, Site specific control of crystallographic grain orientation
ip t
through electron beam additive manufacturing, Materials Science and Technology 31 (8) (2015) 931–938. doi:10.1179/1743284714Y.0000000734.
cr
[38] R. R. Dehoff, M. M. Kirka, F. A. List, III, K. A. Unocic, W. J. Sames, Crystallographic texture engineering through novel melt strategies via electron
us
beam melting: Inconel 718, Materials Science and Technology 31 (8) (2015) 939–944. doi:10.1179/1743284714Y.0000000697.
an
[39] D. G. Duffy, Green’s Functions with Applications, Second Edition, Taylor & Francis Group, 2015.
[40] J. D. Hunt, STEADY-STATE COLUMNAR AND EQUIAXED GROWTH
M
OF DENDRITES AND EUTECTIC, Materials Science and Engineering 65 (1) (1984) 75–83. doi:10.1016/0025-5416(84)90201-5.
d
[41] W. Kurz, C. Bezencon, M. Gaumann, Columnar to equiaxed transition
te
in solidification processing, Science and Technology of Advanced Materials 2 (1, SI) (2001) 185–191. doi:10.1016/S1468-6996(01)00047-X.
Ac ce p
[42] M. A. Groeber, M. A. Jackson, Dream.3d: A digital representation environment for the analysis of microstructure in 3d, Integrating Materials and Manufacturing Innovation 3 (1) (2014) 5. doi:10.1186/2193-9772-3-5.
[43] R. Balluffi, S. Allen, W. Carter, Kinetics of Materials, John Wiley & Sons, 2005.
[44] H. Hu, S. A. Argyropoulos, Mathematical modelling of solidification and melting: a review, Modelling and Simulation in Materials Science and Engineering 4 (4) (1996) 371. [45] V. R. Dave, D. A. Hartman, M. J. Cola, In-Process Quality Assurance for Aerospace Welding, Welding Journal 88 (2) (2009) 28–32.
47
Page 47 of 48
[46] M. Uchic, M. Groeber, M. Shah, P. Callahan, A. Shiveley, M. Scott, M. Chapman, J. Spowart, An Automated Multi-Modal Serial Sectioning
ip t
System for Characterization of Grain-Scale Microstructures in Engineering
Materials, in: DeGraef, M and Poulsen, HF and Lewis, A and Simmons, J
cr
and Spanos, G (Ed.), Proceedings of the 1st International Conference on 3D Materials Science, 2012, pp. 195–202, Seven Springs, PA, Jul. 08-12,
Ac ce p
te
d
M
an
us
2012.
48
Page 48 of 48