# Fileset

[Katagiri_2023_CaseStudies_inThermalAnalysis.pdf](https://mdr.nims.go.jp/filesets/d74ee0f4-bb7e-42b9-942b-0d5695dddf05/download)

## Creator

[Jun Katagiri](https://orcid.org/0000-0002-6399-1951), [Masahiro Kusano,](https://orcid.org/0000-0002-5061-0195), [Sukeharu Nomoto](https://orcid.org/0000-0001-6503-4699), [Makoto Watanabe](https://orcid.org/0000-0002-5064-9583)

## Rights

[Creative Commons BY-NC-ND Attribution-NonCommercial-NoDerivs 4.0 International](https://creativecommons.org/licenses/by-nc-nd/4.0/)

## Other metadata

[Influence of recoil pressure, mushy zone flow resistance and reflectivity on melt pool shape in laser powder bed fusion simulation](https://mdr.nims.go.jp/datasets/68e59dfa-5143-4db0-aa37-979073537ab8)

## Fulltext

Influence of recoil pressure, mushy zone flow resistance and reflectivity on melt pool shape in laser powder bed fusion simulationCase Studies in Thermal Engineering 50 (2023) 103477A2(Contents lists available at ScienceDirectCase Studies in Thermal Engineeringjournal homepage: www.elsevier.com/locate/csiteInfluence of recoil pressure, mushy zone flow resistance andreflectivity on melt pool shape in laser powder bed fusion simulationJun Katagiri ∗, Masahiro Kusano, Sukeharu Nomoto, Makoto WatanabeAdditive Manufacturing Group, Research Center for Structural Materials, National Institute for Materials Science, Sengen1-2-1, Tsukuba, 305-0047, Ibaraki, JapanA R T I C L E I N F OKeywords:Laser powder bed fusionMulti-physics CFD simulationRecoil pressureMushy zone flow resistanceFresnel laser reflection modelMelt pool shapeA B S T R A C TMelt pool dimensions (depth: 𝐷 and width: 𝑊 ) are strongly correlated with defects generatedin laser powder bed fusion. The melt pool dimensions have been evaluated by CFD simulationincluding solid–liquid–gas phase change and laser multiple reflection. Such a multi-physicssimulation requires a large number of parameters. Preliminary simulations for parameteridentification require a lot of time and effort. The parameter sensitivity to the melt pooldimensions is important; however, the systematic evaluations have hardly been found in theliterature. In this study, we performed systematic parametric evaluation for three parametersrelated to recoil pressure, flow resistance in solid–liquid mushy zone and Fresnel laser reflection.As a result, the 𝐷 values increased with increasing recoil pressure, but the 𝑊 values did not.The flow resistance force influenced the velocities in the mushy zone but not the melt pooldimensions. The 𝐷 increased with increasing reflectivity, but the 𝑊 did not. The melt pooldimensions varying recoil pressure, flow resistance force and reflectivity were all inconsistentwith the experimental melt pool dimensions. In order to agree with the experimental melt pooldimensions, the problems to be solved were discussed by comparing the previous studies.1. IntroductionLaser powder bed fusion (L-PBF) is one of the most promising additive manufacturing processes. First, a thin layer of metalpowder is applied to a substrate. A laser moves horizontally and then selectively melts and solidifies the powder bed to produce asliced object derived from an arbitrarily shaped three-dimensional object. A three-dimensional part can be produced by repeatingthe unit process. An advantage of the L-PBF process is that it can produce complex-shaped objects that cannot be produced byconventional metal processing. However, several types of defects occur during L-PBF fabrication, depending on the laser power andlaser scanning speed (the laser processing parameters) [1,2].To identify the laser process parameter that causes the defects, a number of single-track experiments have been carried out,e.g. [2,3]. Additionally, multi-physics Computational Fluid Dynamics (CFD) simulations have been used in previous studies toevaluate the melt pool shape. Such simulations implement multiple physical phenomena such as heat conduction, fluid dynamics,solid–liquid–gas phase change are coupled. The multi-physics simulation has been successfully modelled keyhole-type defect [4],lack-of-fusion type defect during multi-track simulation [5], and balling defect [6]. Some previous studies have reported that the meltpool dimensions obtained from the simulation are in good agreement with those obtained from the experimental measurement [7,8].The multi-physics simulation of the L-PBF process requires various physical models such as recoil pressure due to evaporation [8–10], flow resistance force in solid–liquid coexisting (mushy) zone [11], and Fresnel laser reflection model [12]. Such models include∗ Corresponding author.E-mail address: KATAGIRI.Jun@nims.go.jp (J. Katagiri).vailable online 11 September 2023214-157X/© 2023 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND licensehttp://creativecommons.org/licenses/by-nc-nd/4.0/).https://doi.org/10.1016/j.csite.2023.103477Received 28 April 2023; Received in revised form 31 August 2023; Accepted 7 September 2023https://www.elsevier.com/locate/csitehttp://www.elsevier.com/locate/csitemailto:KATAGIRI.Jun@nims.go.jphttps://doi.org/10.1016/j.csite.2023.103477http://crossmark.crossref.org/dialog/?doi=10.1016/j.csite.2023.103477&domain=pdfhttps://doi.org/10.1016/j.csite.2023.103477http://creativecommons.org/licenses/by-nc-nd/4.0/Case Studies in Thermal Engineering 50 (2023) 103477J. Katagiri et al.fssprpowo22foparameters to be determined. These parameter values have sometimes been determined empirically. The recoil pressure and the flowresistance force are included in the conservation of momentum equation. Therefore, the values of the parameters affect the velocityof the molten metal; in other words, the parameter values will affect the melt pool shape. In addition, the previous experimentalstudies have shown that the laser process parameters extremely influence the melt pool shape [2], which implies an importance ofthe other laser properties such as the laser reflectivity. In this study, the laser multiple reflection is simulated using the ray-tracingmethod. The laser reflectivity in ray-tracing is calculated using the Fresnel model proposed by Cho and Na [12]. The model alsoincludes an empirical parameter. Li et al. have evaluated the effect of the flow resistance force parameter in the mushy zone usingthe multi-physics simulation [13]. They reported that the mushy zone constant influenced the velocity fields inside the melt pooland the surface shape of the top of the melt pool. Alphonso et al. have performed the multi-physics simulations with and withoutrecoil pressure [14] and reported that the melt pool depth obtained from the simulation with recoil pressure was greater than thatwithout recoil pressure. Despite the importance of the three parameters for the melt pool shape, i.e., the defects to be formed, thereis little systematic evaluation of the parameters in the literature. The users of the L-PBF machine should choose the material to beused, the powder properties, the part design, the process parameters and so on. To find the optimum processing parameters forL-PBF, a series of preliminary experiments varying the process parameters requires a lot of time and money. Computer simulationshould be an extremely powerful tool; such a virtual testing tool enables us to efficiently optimise the processing parameters ofL-PBF.In this study, three case studies were carried out to evaluate the effect of recoil pressure, mushy zone flow resistanceorce and laser reflectivity on the melt pool dimensions. In addition, the single-track experiments were also carried out. Theimulation methodology, including its governing equations and the physical models used, and the single-track experiment are brieflyummarised in Section 2. In Section 3, a series of simulations were carried out to vary the three parameters. We found that the recoilressure and the reflectivity parameters affected the melt pool depth, while the flow resistance force parameter did not. The Fresneleflection model parameter was the most influential parameter on the melt pool depth among the three parameters. All the threearameters hardly influenced the melt pool width. Although the three parameters varied significantly, the melt pool dimensionsbtained from the simulations did not agree with those obtained from the single-track experiments. The reason for this discrepancyas then discussed. Finally, the conclusions of this study and the problems to be solved for simulating the melt pool dimensionsbtained from the experimental measurement are described in Section 4.. Materials and methods.1. Discrete element modelling of powder bed formationThe powder bed formed on the substrate is modelled using the Discrete Element Method (DEM) [15]. The DEM is used to simulaterictional rigid particle assemblies. We used the LIGGGHTS (version 3.8.0) to generate the powder bed. More than 80% by volumef the powders used in single-track experiment as described below were between 20 μm and 40 μm in particle diameter, with anaverage diameter of 28.9 μm. To model the particle size distribution, the spherical particles, whose diameter ranged between 20 μmand 40 μm and average diameter was 27.5 μm, were used in the DEM simulation.The Hertz contact model was used in the simulation and its stiffness (Young’s modulus) and Poisson’s ratio were set to 1.3 × 1011N/m and 0.33, respectively. The particle damping coefficient was set to be 0.5 of the particle restitution coefficient. The frictioncoefficients of particle–particle contact and particle–wall contact were both 0.1. The boundary conditions for the six side faces wereset to the f boundary style implemented in LIGGGHTS. The f style boundary basically acts as a fixed smoothed wall boundary, butthe particles are lost if the particle position is outside the domain.The overview of the powder bed formation simulation is as follows. First, the particles were generated to a rectangular areaabove the powder bed to be formed. The particles were deposited under gravity on the powder bed stage, whose width and lengthare 1.0 mm and 0.5 mm respectively. An initial height of the powder bed stage is 0.03 mm. After the gravitational deposition, arigid rectangular wall moves in the 𝑦 direction for 1.0 mm to remove the particles whose height is higher than 0.03 mm.To incorporate the formed powder bed into the multi-physics CFD simulation, a stereolithography (STL) was generated usingthe positions and radii of all particles obtained from the DEM simulation. The powder particles are therefore fixed, which meansthat the particle packing structure does not change during the simulation described below.2.2. Governing equations and physical modelsThe simulation solvers are based on a general purpose CFD software of OpenFOAM (version 8) and partly customised formodelling the L-PBF: TBinterFoam solver (TERRABYTE Co., Ltd., Tokyo, Japan). The short summary of the TBinterFoam solver isdescribed in this section.The governing equations of conservation of energy, conservation of mass and conservation of momentum are expressedby Eqs. (1)–(3) respectively.𝜕(𝜌𝐶𝑝𝑇 )𝜕𝑡+𝜕(𝜌𝑢𝑗𝐶𝑝𝑇 )𝜕𝑥𝑗= −𝜕𝑞𝑗𝜕𝑥𝑗+ ̇𝑞𝑚 + ̇𝑞𝑣, (1)𝜕𝑢𝑖 = 0, (2)2𝜕𝑥𝑖Case Studies in Thermal Engineering 50 (2023) 103477J. Katagiri et al.wcmA(v𝜕(𝜌𝑢𝑖)𝜕𝑡+𝜕(𝜌𝑢𝑖𝑢𝑗 )𝜕𝑥𝑖= −𝑔𝑗(𝑥𝑗 − 𝑟𝑗) 𝜕𝜌𝜕𝑥𝑖+ 𝜕𝜕𝑥𝑗{𝜇(𝜕𝑢𝑖𝜕𝑥𝑗+𝜕𝑢𝑗𝜕𝑥𝑖)}+ 𝐹𝐶𝑆𝐹 ,𝑖 + 𝐹𝑟,𝑖 + 𝑆𝑖, (3)here 𝜌, 𝐶𝑝, 𝑢𝑖, 𝑟𝑗 , 𝑞𝑗 , 𝜇, ̇𝑞𝑚 and ̇𝑞𝑣 are the density, specific heat capacity, flow velocity, reference point, heat flux, viscosity, phasehange heat dissipation term and phase change heat absorption term, respectively. Eqs. (1)–(3) were solved by the finite volumeethod (FVM) in OpenFOAM. FVM is one of the most widely used methods for computational fluid dynamics (CFD) simulation.computational domain is partitioned by a set of polyhedral volume elements called control volumes. The governing equationspartial differential equations) are discretised to algebraic equations. The discretised equations are conservative for each contrololume. The FVM and OpenFOAM software are described in detail in the references cited [16–19].The term 𝐹CSF,𝑖 is the surface tension term using continuum surface force (CSF) model which is expressed as follows.𝐹𝐶𝑆𝐹 ,𝑖 = 𝜎𝜅 𝜕𝐹𝜕𝑥𝑖, (4)where 𝜎 and 𝜅 are the surface tension coefficient and the curvature of the gas–liquid interface, respectively. Note that the interfaceof solid or metal (phase 1) and gas (phase 2) is captured by the function 𝐹 using the Volume-of-Fluid (VOF) method , which isexpressed as follows:𝜕𝐹𝜕𝑡+𝜕𝐹𝑢𝑗𝜕𝑥𝑗= 0. (5)The 𝐹 is the value of the VOF function and ranges between 0 and 1. The gas–liquid interface is represented by the volume ratio ofphase-1 and phase-2 and its dynamics is calculated by solving the advection equation of 𝐹 .The term 𝐹r,𝑖 in Eq. (3) is recoil pressure and is expressed by Eq. (6).𝐹𝑟,𝑖 = −𝛼𝑃0 exp((𝐿𝑣𝑇 − 𝑇𝑣𝑅𝑇𝑇𝑣))𝐧, (6)where, 𝛼, 𝑃0, 𝐿𝑣, 𝑅, 𝑇𝑣, and 𝐧 are an arbitrary constant value, the gas phase pressure, the latent heat of evaporation, the gas constant,the evaporation temperature, and a unit normal vector of the gas-liquid interface.The term 𝑆𝑖 in Eq. (3) represents flow resistance force mushy zone, i.e., solid–liquid coexisting zone and is known asVoller–Prakash model [11]. The Voller–Prakash model is described by Eq. (7).𝑆𝑖 = −𝐶(1 − 𝑔𝑙)2𝑔𝑙3𝑢𝑖, (7)where 𝐶 and 𝑔𝑙 are an arbitrary constant and the volume ratio of the liquid phase, respectively.An enthalpy is calculated using Eq. (8).𝐻 = 𝑔𝑠 ∫𝑇𝑇𝑟𝑒𝑓𝜌𝑠𝐶𝑠𝑑𝑇 + 𝑔𝑙 ∫𝑇𝑇𝑟𝑒𝑓𝜌𝑙𝐶𝑙𝑑𝑇 + 𝜌𝑙𝐶𝑙𝐿, (8)where, 𝑔𝑠, 𝐶𝑠, 𝐶𝑙, 𝑇𝑟𝑒𝑓 , and 𝐿 are the solid fraction, i.e., 𝑔𝑙+𝑔𝑠=1, specific heat of solid phase, specific heat of liquid phase,reference temperature, and latent heat of fusion, respectively. Assuming that the latent heat and the density are independent ofthe temperature, Eq. (9) is derived.𝐻 = 𝑔𝑠𝜌𝑠𝐶𝑠𝑇 + 𝑔𝑙𝜌𝑙𝐶𝑙𝑇 + 𝑔𝑙𝜌𝑙𝐿. (9)The term ̇𝑞𝑚 in Eq. (1) is expressed by following equation.̇𝑞𝑚 = −𝜌𝑙𝐿𝜕𝑔𝑙𝜕𝑡. (10)The enthalpy considering the evaporation is expressed by Eq. (11).𝐻 = 𝑔𝑠 ∫𝑇𝑇𝑟𝑒𝑓𝜌𝑠𝐶𝑠𝑑𝑇 + 𝑔𝑙 ∫𝑇𝑇𝑟𝑒𝑓𝜌𝑙𝐶𝑙𝑑𝑇 + 𝜌𝑙𝐶𝑙𝐿 + 𝑔𝑣 ∫𝑇𝑇𝑟𝑒𝑓𝜌𝑣𝐶𝑣𝑑𝑇 + 𝜌𝑣𝐶𝑣𝐿𝑣, (11)where, 𝑔𝑣, 𝐶𝑣, 𝜌𝑣, and 𝐿𝑣 are the gas fraction, specific heat of evaporation, gas density, and latent heat of evaporation, respectively.Assuming that the latent heat and the density are independent of the temperature, Eq. (12) is derived.𝐻 = 𝑔𝑠𝜌𝑠𝐶𝑠𝑇 + 𝑔𝑙𝜌𝑙𝐶𝑙𝑇 + 𝑔𝑣𝜌𝑣𝐶𝑣𝑇 + 𝑔𝑙𝜌𝑙𝐿 + 𝑔𝑣𝜌𝑣𝐿𝑣. (12)The term ̇𝑞𝑣 in Eq. (1) is expressed by following equation.̇𝑞𝑣 = −(𝜌𝑣𝐶𝑣(𝑇 − 𝑇𝑟𝑒𝑓 ) + 𝜌𝑣𝐿𝑣) 𝜕𝑔𝑣𝜕𝑡. (13)Heat flux (𝑞) between the metal surface and the gas phase is expressed by Eq. (14).𝑞 = −ℎ(𝑇 − 𝑇𝑔) − 𝜖ℎ𝜎𝑏(𝑇 4 − 𝑇𝑔4), (14)where ℎ, 𝑇𝑔 , 𝜖ℎ and 𝜎𝑏 are the heat transfer coefficient, gas phase temperature, surface emissivity and Stefan–Boltzmann constant(5.67×10−8 Wm−2K−4) respectively.3Case Studies in Thermal Engineering 50 (2023) 103477J. Katagiri et al.Fig. 1. A schematic of Fresnel reflection model.Various moving heat source models have been used in the simulation of the L-PBF process, such as the application of a Gaussianheat source to the surface [8], the volumetric heat source model in which the laser energy is absorbed depending on the depthfrom the metal surface [20,21], and the ray tracing method [5,6]. The ray tracing method is one of the solvers used in this study.The laser beam is assumed to consist of a bundle of laser beams. The ray tracing method tracks each ray. Fig. 1 shows a schematicof the ray tracing method. When the laser reaches the metal surface, thermal energy is applied to the metal surface. The thermalenergy applied has been calculated in inverse proportion to the laser reflectivity of the metal surface. The laser beam repeats themirror reflections until (1) the applied energy became extremely small or (2) the laser beam did not reach the metal surface. Thereflectivity (𝑅𝐿) is calculated using the Fresnel reflection model in the ray tracing method [12].𝑅𝐿 = 12(1 + (1 − 𝜖 cos 𝜃)21 + (1 + 𝜖 cos 𝜃)2+ 𝜖2 − 2𝜖 cos 𝜃 + 2 cos2 𝜃𝜖2 + 2𝜖 cos 𝜃 + 2cos2𝜃), (15)where, 𝜃 and 𝜖 are the laser incident angle and an arbitrary constant, respectively.The laser energy distribution implemented in the TBinterFoam solver is a Gaussian distribution. The heat flux (𝑞𝑏) which isexpressed as follows.𝑞𝑏 = 𝐶(𝑃 ) exp(−𝑟2∕𝑟𝐿22𝜎𝐺2), (16)where, 𝐶(𝑃 ), 𝑟𝐿, and 𝜎𝐺 are a scaling factor depending on the laser power (𝑃 ), the laser spot radius (40×10−6 m in this study),and the standard deviation of the Gaussian distribution (1.0 in this study). Note that the formula for 𝐶(𝑃 ) is not described in thevendor’s manual.The explicit Euler time integration scheme, the linear interpolation scheme for the Laplacian and gradient discretisation, and thesecond-order linear upwind and total variation diminishing (TVD) schemes for the divergent discretisation were used. Note that thethird-order linear upwind scheme is often used in CFD simulation. We chose the TVD scheme because the preliminary simulationusing the TVD scheme was more stable compared to the third-order linear upwind scheme.The Gauss–Seidel method was used to calculate the flow and temperature fields. The tolerance value of the flow velocityand temperature was 1.0×10−6. The incomplete Cholesky preconditioned conjugate gradient method was used for the pressurecalculation. Considering the balance between calculation accuracy and speed, the residual value for the iteration was 1.0×10−6 toensure the stability of the calculation.Note that the governing equations and physical models used are the same as those used in the previous studies of L-PBFsimulation [3–8,14,20–23]. The dynamics of the gas phase is sometimes not considered in the previous studies; on the other hand,both the dynamics of the fluid and the gas phase are considered in this study.2.3. Physical properties and calculation parameters usedThe material used in this study is Inconel738LC which is a nickel based alloy. The material properties and simulation parametersused are listed in Table 1. Note that the material properties are based on the calculation of CALPHAD software: JMatPro(www.sentesoftware.co.uk/jmatpro (accessed 7 April 2023)).The dynamics of the melt pool is calculated using Eq. (3). Eq. (3) includes the source terms of Eqs. (4), (6) and (7), which meansthat the parameters in these equations can influence the melt pool dimensions. Eq. (15) is not included in Eq. (3); however, Eq. (15)directly influences the heat flux applied to the metal. Hence, we have performed a series of parameter sensitivity analyses for thethree parameters.4https://www.sentesoftware.co.uk/jmatproCase Studies in Thermal Engineering 50 (2023) 103477J. Katagiri et al.ospmwpaocitdoteadmpfbmsTCtTable 1Materials properties of Inconel738LC and the simulation parameters used in this study.Property or parameter Value [Unit]Density of solid 8820 [kg/m3]Density of liquid 8220 [kg/m3]Density of Gas 1.176 [kg/m3]Solidus temperature 1371.15 [K]Liquidus temperature 1620.15 [K]Evaporation temperature 3034.0 [K]Thermal conductivity 46.1 [W/(m K)]Latent heat of fusion 25,000 [J/kg]Latent heat of evaporation 734,000 [J/kg]Specific heat capacity of metal 710 [J/kg K]Specific heat capacity of gas 1007 [J/kg K]Surface tension coefficient (𝜎) 1.2 [N/m]Liquid viscosity (𝜇) 0.005 [Pa s]Gas viscosity 1.865 × 10−5 [Pa s]Liquid–gas heat transfer coefficient 40 [W/(m2 K)]Emissivity of thermal radiation 0.3 [–]Laser power (𝑃 ) 300 [W]Laser scanning speed (𝑉 ) 1.0 [m/s]Laser diameter 80.0 × 10−6 [m]Laser type Parallel beam ray [–]Laser energy distribution Gaussian [–]Standard deviation of Gaussian distribution 1.0 [–]Value of 𝛼 in Eq. (6) 0.05 - 500 [–]Value of 𝐶 in Eq. (7) 10−1 - 109 [kg/(m3 s)]Value of 𝜖 in Eq. (15) 0.2 - 1.0 [–]2.4. SimulationFig. 2(a) shows the snapshot at the start of the simulation. We used a rectangular simulation model with width, depth and heightf 0.5 mm (−0.25 mm ≤ 𝑥 ≤ 0.25 mm), 1 mm (0 mm ≤ 𝑦 ≤ 1.0 mm) and 0.3 mm (−0.2 mm ≤ 𝑧 ≤ 0.1 mm) respectively. The gridize used must affect the simulation results. We have performed preliminary simulations varying the grid size; as a result, the meltool dimensions obtained from the simulation became almost the same when the grid size was smaller than 5 μm. Hence, hexagonalesh with the grid size of 5 μm was used in this study. Note that the mesh size was changed to 2.5 μm for a rectangular regionithin the simulation model of 0.4 mm width, 0.8 mm depth and 0.15 mm height. This resulted in a total of 5760,000 cells. Theast studies have used the same or similar grid size; e.g., 3 μm [5], 4 μm [4], 4.5 μm [7], and 5 μm [14].As a typical room temperature, the initial temperatures of metal and atmosphere were set at 298.15 K. For the thermal simulation,constant temperature of 298.15 K was applied to the top 𝑥–𝑦 plane, and an adiabatic boundary condition was applied to thether five planes. For the pressure, an atmospheric pressure (101,325 Pa) was adopted for the upper 𝑥–𝑦 plane, and a zero-gradientondition was adopted for the other five sides. For the flow velocity, the pressureInletOutletVelocity boundary condition implementedn OpenFOAM was adopted for the upper 𝑥–𝑦 plane. The pressureInletOutletVelocity condition assumed a zero gradient condition whenhe flow direction is outward. When the flow direction is inward, the velocity is obtained from the flow with the specified inletirection. The zero-gradient condition was applied to the other five sides.Figs. 2 (b) to (f) show the snapshots during the example case of the simulation for 𝑃 = 300 W and 𝑉 = 0.9 m/s. A starting pointf the laser is indicated by white circle in Fig. 2(a). The laser irradiation direction was vertical downwards. The laser moved towardshe 𝑦 direction for 0.7 mm. The laser process parameters were set to 𝑃 = 300 W and 𝑉 = 1.0 m/s to compare with the single-trackxperiment as described below. When 𝑉 = 1.0 m/s, the simulation time was 0.0007 s. The computation time for a typical case wasbout 60 h on a 96-core parallel computation using an Intel Xeon Platinum 8268 processor (clock speed: 2.9 GHz).After the simulation, the temperature fields were visualised in the cross section of the centre of the model: 0.5 mm in the 𝑦irection. Isolines of solidus and liquidus temperatures were plotted on the visualised image. The dimensions of the melt pool wereeasured from the visualised image. The melt pool depth obtained from the simulation is defined as the length between the deepestoint of the liquidus temperature isoline and the top of the substrate. Fig. 3(a) shows a schematic of the melt pool measurementor the CFD simulation and the single track experiment. The melt pool width obtained from the simulation is defined as the lengthetween the minimum and maximum points of the liquidus temperature along the 𝑥 direction. The melt pool depth and width wereeasured in the similar way using the melt pool visually recognised from the cross-sectional image.The melt pool obtained from the single-track experiment forms a semicircle on the substrate, whereas that obtained from theimulation does not. The semicircle forms because the surface tension of the melt attracts the surrounding unmelted powder particles.he dynamics of the solid particles is not coupled in the present simulation; therefore the semicircle shape cannot be formed in theFD simulation. In addition, the time history of the temperature at a measurement point shown in Fig. 3(b) was calculated from5he simulation data.Case Studies in Thermal Engineering 50 (2023) 103477J. Katagiri et al.Fig. 2. Snapshots during a sample simulation for 𝑃 = 300 W and 𝑉 = 0.9 m/s at 𝑡 = 0 s (a), 0.00001 s (b), 0.0001 s (c), 0.00025 s (d), 0.0005 s (e), and0.0007 s (f).2.5. Single-track experimentA procedure for the single-track experiment has been described in detail in the cited Refs. [24]. A brief summary of theexperimental procedure is described here.The powder and substrate material used in this study was Inconel738LC. The average particle diameter was 28.9 μm. The scanningtrack was cut in the centre of the track line with a high precision cutting machine. Resin was embedded in the cross section andpolished using a polishing machine. Note that a colloidal silica solution was added during the final polishing process. The crosssection was then washed with water and allowed to dry sufficiently. The surface of the cross section was examined by scanningelectron microscopy (SEM). The width and depth the melt pool was measured from the SEM image as shown in Fig. 3(a). Thistechnique is widely used in the relevant research community to quantify the melt pool dimensions.3. Results and discussions3.1. Case study 1: Influence of recoil pressureFig. 4 shows the time evolution of the temperature for the simulations with 𝛼 = 0.05, 0.5, 5, 50 and 500 respectively. Thetemperature becomes high as the laser source approaches the measurement point and then decreases as the laser source recedes.Such temperature trends have been observed in previous studies [5,7,22]. A maximum temperature occurred at about 𝑡 = 0.0003 sand becomes low as a function of 𝛼, except for 𝛼 = 0.05. The recoil pressure is the force pushing down the metal surface caused bythe ejection of gas particles due to the evaporation. The metal surface is excavated by the recoil pressure. The molten metal massarose from the excavation flows in the opposite direction to the laser movement and accumulates backwards.6Case Studies in Thermal Engineering 50 (2023) 103477J. Katagiri et al.Fig. 3. A schematic of the measurement of melt pool dimensions (a) and the temperature (b).Fig. 4. Time evolution of the temperature for the simulations varying the parameter 𝛼 in recoil pressure arising when the metal vaporisation (Eq. (6)).Fig. 5 shows the metal surface profiles (isoline of VOF function value of 0.5) at 𝑡 = 0.0006 s for 𝛼 = 0.05, 0.54, 5, 50, and 500.The excavated depth at about 𝑦 = 0.75 mm is greater as 𝛼 increases; the trend is similar to the previous study [14]. According to thestudy [14], the melt pool depth obtained from the simulation with recoil pressure is greater than that without recoil pressure. When𝛼 was 0.05, the metal mass hardly accumulated backwards; in other words, the powder bed and the upper part of the substrateappeared to be scraped off. The evaporated volume (mass) is converted to the recoil pressure in the CFD simulation. When the 𝛼 issmall, the recoil pressure is underestimated. As a result, the powder bed and the top of the substrate appeared to disappear.The 𝛼 = 0.54 is a commonly used value in the past studies and is determined by considering the Knudsen layer duringevaporation [9,10,25]. In the CFD simulation algorithm, the recoil pressure 𝐹r is calculated on the basis of 𝑚̇𝑣 ⋅ 𝑣𝑇 , where 𝑚̇𝑣and ̇𝑣𝑇 are the evaporated mass and the velocity at the end of the Knudsen layer [26]. Assuming that 𝑣𝑇 is equal to the local soundvelocity [9,10], Anismov [9] obtained 𝛼 = 0.54. The molten mass moved backwards and the surface height returned to the top of thesubstrate. When 𝛼 ≥ 5, the melt accumulate+d backwards at 1.5 × 10−4 m < 𝑦 < 3.0 × 10−4 m depending on 𝛼. However, the shapeof the melt pool above the substrate often forms a semicircle or a segmented circle whose area is larger than the semicircle in the7Case Studies in Thermal Engineering 50 (2023) 103477J. Katagiri et al.Fig. 5. The isolines of VOF function value of 0.5 at 𝑡 = 0.0006 s for the recoil pressure parameter: 𝛼 = 0.05 (black), 0.54 (red), 5 (light-green), 50 (blue), and500 (light-blue). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)Fig. 6. Relationship between recoil pressure parameter: 𝛼 in Eq. (6) and the melt pool dimension: melt pool depth 𝐷 and width 𝑊 obtained from the simulationand the single-track experiment under 𝑃 = 300 W and 𝑉 = 1.0 m/s.single-track experiment, as shown in Fig. 3(a). The semicircle did not form in the CFD simulations because the semicircle should beformed by two mechanisms: (1) the backward accumulation and (2) the attraction of powder particles due to surface tension. Duringthe particle melting process, the molten liquid should attract and absorb the (unmelted) powder particle by the surface tension. Incontrast, the powder particles in the area near the scanned track are scattered due to the recoil pressure induced by evaporation:the spattering phenomenon [27]. The size and shape of the semicircle on the substrate should depend on the balance between theattracted and absorbed particles and the scattered particles. To reproduce the semicircular shape in the CFD simulation, the solidparticle dynamics should be coupled with the CFD simulation, e.g. DEM-CFD coupling [28,29]. Note that the recent multi-physicssimulations using DEM-CFD coupling [30,31] and smoothed particle hydrodynamics modelling [32] have incorporated the motionof the solid particle assembly into the multi-physics simulation of L-PBF.Fig. 6 shows the relationship between the 𝛼 and the melt pool dimension: the melt pool depth 𝐷 and width 𝑊 obtained from theCFD simulation and the single-track experiment. As described above, 𝐷 increases with increasing 𝛼 due to the surface excavationcaused by the recoil pressure. On the other hand, 𝑊 is around the laser spot diameter (80 μm) regardless of 𝛼. Since the laserbeam direction is perpendicular downward, the surface planes near the edges of the excavated region are parallel to the laser beamdirection. In such a case, the surface plane hardly absorbs the thermal energy from the laser irradiation according to the Fresnelmodel; as a result, the measured 𝑊 values were almost equal to the laser spot diameter. However, the 𝐷 and 𝑊 values obtained fromthe simulation are different from those obtained from the single-track experiments, which indicates that the melt pool dimensionin the single-track experiment can hardly be simulated just by varying 𝛼.3.2. Case study 2: Influence of mushy zone flow resistanceFig. 7(a) shows the time evolution of the temperature for the simulations with 𝐶 = 10−1, 101, 105, and 107 respectively. Asshown in the figure, the general trend of temperatures for different 𝐶 is almost the same. The 𝐶 is the influencing parameter of thesolid–liquid coexistence region, which means that the 𝐶 is valid between the solidus and liquidus temperature. The temperaturerange is known as brittleness temperature range (BTR). The temperature falls into the BTR at three times: 𝑡 = 0.00028 s (P1),0.00035 s (P2) and 0.00055 s ≤ 𝑡 ≤ 0.00062 s (P3). The temperatures at P1, P2 and P3 are in the solid–liquid translation phase,in the gas phase due to evaporation and in the liquid–solid solidification phase respectively. Eq. (7) is not valid for P2 because8Case Studies in Thermal Engineering 50 (2023) 103477J. Katagiri et al.wwfwFig. 7. Time evolution of the temperature for the simulation varying the mushy zone flow resistance parameter in Eq. (7): 𝐶 = 10−1, 101, 105, 107, and 109.Eq. (7) is for the solid–liquid mushy zone. In the P1 state, the temperature rises rapidly for a short time. As the melting process inP1 stops for the short time, 𝐶 has little effect on the temperature history. The solid–liquid mushy zone at P3 exists for a longer timecompared to the P1 state.Fig. 8 shows the velocity fields at 𝑡 = 0.0006 s in the 𝑦–𝑧 cross section for 𝐶 = 10−1, 101, 105, 107 and 109. The colour bar inthe figure represents the cell velocity component of the 𝑦 direction. The solid white lines in the figures show the isolines of thesolidus and liquidus temperatures. The area between the two white lines is the mushy zone, so the value of 𝐶 affects the velocityin this area. The molten metal mass accumulates backwards and therefore moves in the negative direction 𝑦. Note that the positivevelocities, i.e. the red coloured region, occur in the mushy zone because the accumulated mass forms a convective flow. As shownin the figure, the flow velocities in the mushy zone become small as 𝐶 increases.Fig. 9 shows the effect of 𝐶 on the melt pool dimensions. As shown in the figure, 𝐶 has little effect on the depth and width ofthe melt pool. Voller and Prakash, who developed Eq. (7), proposed the reference values of 𝐶: 1.6×103 and 1.6×104 [11]. In theprevious studies of the L-PBF simulations, the 𝐶 values were 5.57 × 106 [8], 1014 [6] and 105 − 107 [5].Note that the value of 𝐶 cannot be varied in the experiment. The 𝐶 depends on 𝑔l; however, 𝑔l: the volume of the fully moltenmetal region is extremely difficult to measure in the experiment. According to the previous study, the 𝐶 value influences the surfaceshape and the internal velocity field of the melt pool [13]. In the present study, 𝐶 varied widely between 10−1 and 109; however,the melt pool dimensions were almost constant regardless of 𝐶.As discussed in [11], an applicability of the Kozeny–Carman equation should be evaluated. The correlation between 𝐶 and themicro-structural parameters originally included in the Kozeny–Carman equation should be clarified. Voller and Prakash focusedon a pressure drop due to porosity (i.e. liquid fraction) in the Kozeny–Carman equation and derived Eq. (7). They changed thepressure–porosity relation to a force–porosity relation by the parameter 𝐶. The Kozeny–Carman equation includes the micro-structural parameters such as the specific surface area, the hydraulic tortuosity (𝜏ℎ) which characterises the streamline length of porethroats in porous media, and the shape factor which characterises the shape of the cross-section of the porous media perpendicularto the flow direction [33]. The reference values of 𝐶 can be estimated from the micro-structural parameters, but the values ofthe micro-structural parameters depend on other factors such as the physical phenomena being simulated. A universal referencevalue of 𝐶 is difficult to determine, so the 𝐶 values have to be determined by trial and error. The reference value of 𝐶 consideringthe L-PBF process will be determined in our future study. Furthermore, the applicability of the Kozeny–Carman equation should beconsidered, as the Kozeny–Carman equation is used for the Darcy flow which is an extremely low Reynolds number flow. During theL-PBF process, the molten metal mass may become partially inertial flow. In the inertial flow regime, the porous flow is treated as aForchheimer flow [34]. The Ergun equation [35] is a commonly used empirical equation to model the pressure–porosity relationshipfor Forccherimer flow. The flow resistance force using the Ergun equation should be included in the simulation and compared withthat using the Kozeny–Carman equation.3.3. Case study 3: Influence of laser reflectivityFig. 10 shows the relationship between a laser incident angle and the reflectivity values calculated by Eq. (15) for differentvalues of 𝜖. The reflectivity value becomes small as the 𝜖 value increases.Fig. 11 shows the temperature fields in the 𝑥–𝑧 cross section at the centre of the scan track (𝑦 = 0.5 mm) for the simulationith 𝜖 = 0.2, 0.6 and 1.0. The solid white lines are the solidus and liquidus temperatures respectively. The reflectivity decreasesith increasing 𝜖. The metal surface absorbs more thermal energy for larger 𝜖; as a result, the melt pool depth becomes large as aunction of 𝜖.The effect of 𝜖 on the melt pool dimensions is shown in Fig. 12. It is quantitatively clarified that the melt pool depth increasesith increasing 𝜖; however, the 𝐷 obtained from the simulation differs from that obtained from the single-track experiment.9Case Studies in Thermal Engineering 50 (2023) 103477J. Katagiri et al.Fig. 8. Fields of the cell velocity component in 𝑦 direction at 𝑡 = 0.0006 s in the 𝑦–𝑧 cross-section for the mushy zone flow resistance parameter: 𝐶 = 10−1,101, 105, 107, and 109.Fig. 9. Relationship between the mushy zone flow resistance parameter:𝐶 and the melt pool dimensions: melt pool depth 𝐷 and width 𝑊 obtained from thesimulation and the single-track experiment under 𝑃 = 300 W and 𝑉 = 1.0 m/s.The Fresnel model calculates the reflectivity by geometrical reflection analysis, which means that the simulation using the ray-tracing with the Fresnel model is difficult to reproduce the 𝐷 in the single-track experiment. Conventional range of 𝜖 values isbetween 0.2 and 0.25. The direction of the laser irradiation is downward in the 𝑧 direction. The 𝜃 was 0 degrees when the laser wasirradiated in the horizontal plane; in this case the reflectivity is about 0.67 and 0.61 for 𝜖 of 0.2 and 0.25 respectively. Accordingto the previous studies, the melt pool dimensions obtained from the simulations using a constant reflectivity value of 0.7 [23] and10Case Studies in Thermal Engineering 50 (2023) 103477J. Katagiri et al.Fig. 10. Relationship between a laser incident angle and the reflectivity values calculated by Fresnel laser reflection model:Eq. (15) for different values of anempirical parameter: 𝜖.Fig. 11. Temperature fields in the 𝑥–𝑧 cross section at the centre of the scan track (𝑦 = 0.5 mm) for the simulation with the parameter in Fresnel reflectionmodel (Eq. (15)): 𝜖 = 0.2, 0.6 and 1.0. The solid white lines are the solidus and liquidus temperatures respectively.Fig. 12. Relationship between the parameter in Fresnel reflection model (Eq. (15): 𝜖 and the melt pool dimensions: melt pool depth 𝐷 and width 𝑊 obtainedfrom the simulation and the single-track experiment at 𝑃 = 300 W and 𝑉 = 1.0 m/s.11Case Studies in Thermal Engineering 50 (2023) 103477J. Katagiri et al.c𝜎rcrtwrsptsudst0.8 [4] agreed well with those obtained from the experiments. The melt pool dimensions obtained in this study are significantlydifferent from those obtained in the single-track experiment. The physical properties of Inconel738LC in this study are based onCALPHAD and are similar to the previous study for the same material [3], which means that the other parameters are importantto simulate the melt pool dimensions. One possibility is the shape of the Gaussian distribution. The following form of Gaussiandistribution has been used in the field of laser technology and has been used in the L-PBF simulations in numerous studies in thepast, e.g. [3,5–8].𝑞𝑏 =2𝑃𝐴𝜋𝑟𝐿2exp(−2 𝑟2𝑟𝐿2), (17)where 𝐴 is the absorptivity. Since the transmissivity of metal is zero, 𝐴 = 1 - 𝜏. The 𝑞𝑏 is linearly correlated with 𝐴 in Eq. (17).The thermal energy applied to the metal surface varies as a function of 𝐴, which implies that the laser reflectivity must influencethe shape of the melt pool. Eq. (16) is equivalent to Eq. (17) if 𝜎 = 0.5. The 𝜎 in this study was 1.0. The laser energy is moreoncentrated when the 𝜎 is smaller, which implies that the melt pool must be deeper when 𝜎 < 1.0. The simulations with differentvalues will be performed in our future study.On the other hand, Ren et al. have reported that the absorptivity obtained from the simulation using temperature-dependenteflectivity model, so-called Hagan–Ruben model was 0.404 (reflectivity: 0.596) and the effective absorptivity obtained from thealorimetric measurements was 0.424 (reflectivity: 0.576) under the condition of 𝑃 = 430 W and 𝑉 = 0.6 m/s. They have alsoeported that the time-averaged absorptivity was 0.094 (reflectivity: 0.906) for the simulation using the Fresnel reflection model,he form of which is slightly different from Eq. (15). According to Ren et al. [6], the Hagan-Ruben model is reasonable comparedith the Fresnel model. Different reflectivity models and their values have been used in the previous studies and the value of theeflectivity significantly influences the melt pool dimensions obtained from the simulation.The melt pool width: 𝑊 increases slightly with increasing 𝜖. The 𝑊 obtained from the simulation does not agree with thatobtained from the single-track experiment. The discrepancy may be caused by the effect of reflectivity perpendicular to the laserirradiation direction: around zero degrees of the laser incident angle. Fig. 10 shows that the reflectivity value varies significantlywith 𝜖. Since the laser irradiation direction is vertically downwards, the number of laser beams perpendicular to the laser directionshould be extremely small. Although the reflectivity values vary significantly, the melt pool width in the horizontal direction wasdifficult to expand due to the small number of rays perpendicular to the laser direction. However, the melt pool width around theupper part of the substrate is often larger than that inside the melt pool. The larger 𝑊 in the upper part may reflect the differencein melting behaviour between the powder bed and the substrate. To clarify the reason and to simulate the melt pool dimensionsobtained from the single-track experiment, further studies for the simulations using constant reflectivity and the Hagen–Ruben modelshould be required.4. ConclusionsIn this study, a series of multi-physics CFD simulations were carried out to evaluate the influence of the parameters for therecoil pressure, the flow resistance force in the solid–liquid mushy zone and the Fresnel reflection model. The simulation resultswere quantitatively compared with the single-track experiments, and the following conclusions were drawn.• The melt pool depth increased with increasing parameter 𝛼 for the recoil pressure, while the melt pool width did not. Themelt pool dimension obtained from the simulations with varying 𝛼 did not agree with those obtained from the single-trackexperiment.• The evaporated mass is converted to recoil pressure in the CFD algorithm. When the large 𝛼 was used, the metal surface wasdeeply excavated and then the molten mass flowed backwards and accumulated. When the extremely low 𝛼, the molten metalmass hardly accumulated backwards; as a result, the powder bed and the upper part of the substrate disappeared after passingthe laser.• The melt pool depth and width were almost constant irrespective of the parameter 𝐶 for the flow resistance force in the mushyzone, and did not agree with that obtained from the single-track experiment.• The melt pool depth increased with increasing parameter 𝜖 in the Fresnel reflection model, while the melt pool width did not.The melt pool depth and width obtained from the simulations with varying 𝜖 did not agree with that obtained from the singletrack experiment.The material used in this study is Inconel738LC; however, the melting and solidification behaviour of the metallic materialshould be similar. Therefore, the results of this study will be applicable to the L-PBF for the metallic materials.To produce a defect-free object using the L-PBF process, engineers should consider the effects of laser processing conditions (laserower, scan speed) and powder characteristics (particle size distribution, particle shape, packing homogeneity). It is quite difficulto optimise such processing parameters by experiment alone. In addition, the L-PBF process is capable of fabricating complex-haped objects compared to the conventional fabrication processes. To take advantage of the features of L-PBF, it is essential tose computer-aided engineering (CAE) technologies, such as optimal design of the lattice structure, evaluation of residual stresses,efects and the microstructure to be formed. The multi-physics simulation is used in this study to evaluate defect formation. Foruch simulations to become a virtual testing tool, the melt pool dimensions should be reproduced by the simulation. In the light ofhe above discussion, the problems to be solved are as follows.12Case Studies in Thermal Engineering 50 (2023) 103477J. Katagiri et al.tetmmCpeFDiDAPJTR• To simulate a semicircular shape often observed in the experiment, not only the parameter tuning but also the solid particledynamics such as coupling CFD and DEM and SPH should be incorporated to the simulation.• The effect of microstructural parameters in the Kozeny–Carman equation, which is the basis of the Voller-Prakash model,should be evaluated. The flow resistance force model in the mushy zone should be extended for the inertial flow regime.• The laser reflection model should include not only simple geometric reflection but also the temperature dependency.• The standard deviation of the Gaussian distribution should be evaluated.In this study, the processing parameter was a single condition: 𝑃 = 300 W and 𝑉 = 1.0 m/s. In order to verify the validity ofhe simulation, the authors have started a series of simulations where 𝑃 and 𝑉 are varied. The simulation results will be publishedlsewhere.Not only the melt pool shape and defect formation, but also the microstructure is extremely important for the quality ofhe manufactured part. Pioneering work on the coupling of CFD with microstructure formation analysis, such as the phase fieldodel [36] and the cellular automaton model [37], has been published in recent years. We will develop the coupling of CFD andicrostructure formation model for microstructure control.RediT authorship contribution statementJun Katagiri: Conceptualization, Methodology, Validation, Formal analysis, Investigation, Data curation, Writing – original draftreparation, Writing – review & editing, visualization. Masahiro Kusano: Software, Validation, Formal analysis, Writing – review &diting. Sukeharu Nomoto: Methodology, Validation, Formal analysis, Writing – review & editing. Makoto Watanabe: Validation,ormal analysis, Writing – review & editing, Supervision, Project administration, Funding acquisition.eclaration of competing interestThe authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, ornterpretation of data; in the writing of the manuscript; or in the decision to publish the results.ata availabilityThe data that has been used is confidentialcknowledgementsFinancial support from the Council for Science, Technology, and Innovation (CSTI), Cross-ministerial Strategic Innovationromotion Programme (SIP), ‘‘Materials Integration for revolutionary design system of structural materials’’ (Funding agency:ST) and Innovative Science and Technology Initiative for Security, Grant Number JPJ004596 implemented by the Acquisition,echnology, and Logistics Agency(ATLA)’’.eferences[1] T. DebRoy, H.L. Wei, J.S. Zuback, T. Mukherjee, J.W. Elmer, J.O. Milewski, A.M. Beese, A. Wilson-Heid, A. De, W. Zhang, Additive manufacturing ofmetallic components – process, structure and properties, Prog. Mater. Sci. 92 (2018) 112–224, http://dx.doi.org/10.1016/j.pmatsci.2017.10.001.[2] R. Seede, D. Shoukr, B. Zhang, A. Whitt, S. Gibbons, P. Flater, A. Elwany, R. Arroyave, I. Karaman, An ultra-high strength martensitic steel fabricatedusing selective laser melting additive manufacturing: Densification, microstructure, and mechanical properties, Acta Mater. 186 (2020) 199–214, http://dx.doi.org/10.1016/j.actamat.2019.12.037.[3] D. Grange, A. Queva, G. Guillemot, M. Bellet, J.D. Bartout, C. Colin, Effect of processing parameters during the laser beam melting of Inconel 738:Comparison between simulated and experimental melt pool shape, J. Mater Process. Technol. 289 (2021) http://dx.doi.org/10.1016/j.jmatprotec.2020.116897.[4] L. Wang, Y. Zhang, H.Y. Chia, W. Yan, Mechanism of keyhole pore formation in metal additive manufacturing, npj Comput. Mater. 8 (2022)http://dx.doi.org/10.1038/s41524-022-00699-6.[5] M. Bayat, S. Mohanty, J.H. Hattel, Multiphysics modelling of lack of fusion voids formation and evolution in In718 made by multitrack/multi-layer L-PBF,Int. J. Heat Mass Transfer 139 (2019) 95–114, http://dx.doi.org/10.1016/j.ijheatmasstransfer.2019.05.003.[6] Z. Ren, D.Z. Zhang, G. Fu, J. Jiang, M. Zhao, High-fidelity modelling of selective laser melting copper alloy: Laser reflection behavior and thermal-fluiddynamics, Mater. Des. 207 (2021) http://dx.doi.org/10.1016/j.matdes.2021.109857.[7] P. Ninpetch, P. Kowitwarangkul, S. Mahathanabodee, P. Chalermkarnnon, P. Rattanadecho, Computational investigation of thermal behavior and moltenmetal flow with moving laser heat source for selective laser melting process, Case Stud. Therm. Eng. 24 (2021) http://dx.doi.org/10.1016/j.csite.2021.100860.[8] W. Wang, W. Lin, R. Yang, Y. Wu, J. Li, Z. Zhang, Z. Zhai, Mesoscopic evolution of molten pool during selective laser melting of superalloy Inconel 738at elevating preheating temperature, Mater. Des. 213 (2022) http://dx.doi.org/10.1016/j.matdes.2021.110355.[9] S.I. Anisimov, Vaporization of metal absorbing laser radiation, J. Exp. Theor. Phys. 27 (1968) 182–183.[10] A. Klassen, T. Scharowsky, C. Körner, Evaporation model for beam based additive manufacturing using free surface lattice Boltzmann methods, J. Phys.D: Appl. Phys. 47 (2014) http://dx.doi.org/10.1088/0022-3727/47/27/275303.[11] V.R. Voller, C. Prakash, A fixed grid numerical modelling methodology for convection–diffusion mushy region phase-change problems, Int. J. Heat MassTransfer 30 (1987) 1709–1719.[12] J.H. Cho, S.J. Na, Implementation of real-time multiple reflection and fresnel absorption of laser beam in keyhole, J. Phys. D: Appl. Phys. 39 (2006)5372–5378, http://dx.doi.org/10.1088/0022-3727/39/24/039.13http://dx.doi.org/10.1016/j.pmatsci.2017.10.001http://dx.doi.org/10.1016/j.actamat.2019.12.037http://dx.doi.org/10.1016/j.actamat.2019.12.037http://dx.doi.org/10.1016/j.actamat.2019.12.037http://dx.doi.org/10.1016/j.jmatprotec.2020.116897http://dx.doi.org/10.1016/j.jmatprotec.2020.116897http://dx.doi.org/10.1016/j.jmatprotec.2020.116897http://dx.doi.org/10.1038/s41524-022-00699-6http://dx.doi.org/10.1016/j.ijheatmasstransfer.2019.05.003http://dx.doi.org/10.1016/j.matdes.2021.109857http://dx.doi.org/10.1016/j.csite.2021.100860http://dx.doi.org/10.1016/j.csite.2021.100860http://dx.doi.org/10.1016/j.csite.2021.100860http://dx.doi.org/10.1016/j.matdes.2021.110355http://refhub.elsevier.com/S2214-157X(23)00783-9/sb9http://dx.doi.org/10.1088/0022-3727/47/27/275303http://refhub.elsevier.com/S2214-157X(23)00783-9/sb11http://refhub.elsevier.com/S2214-157X(23)00783-9/sb11http://refhub.elsevier.com/S2214-157X(23)00783-9/sb11http://dx.doi.org/10.1088/0022-3727/39/24/039Case Studies in Thermal Engineering 50 (2023) 103477J. Katagiri et al.[13] E. Li, Z. Zhou, L. Wang, H. Shen, R. Zou, A. Yu, Particle scale modelling of melt pool dynamics and pore formation in selective laser melting additivemanufacturing, Powder Technol. 397 (2022) 117012, http://dx.doi.org/10.1016/j.powtec.2021.11.056.[14] W. Alphonso, M. Baier, S. Carmignato, J. Hattel, M. Bayat, On the possibility of doing reduced order, thermo-fluid modelling of laser powder bed fusion(L-PBF) – assessment of the importance of recoil pressure and surface tension, J. Manuf. Process. 94 (2023) 564–577, http://dx.doi.org/10.1016/j.jmapro.2023.03.040.[15] P.A. Cundall, O.D.L. Strack, A discrete element model for granular assemblies, Géotechnique 29 (1979) 47–65.[16] H. Versteeg, W. Malalasekera, An Introduction to Computational Fluid Dynamics: The Finite Volume Method, Pearson Education Limited, 2007.[17] H.G. Weller, G. Tabor, H. Jasak, C. Fureby, A tensorial approach to computational continuum mechanics using object-oriented techniques, Comput. Phys.12 (6) (1998) 620–631, http://dx.doi.org/10.1063/1.168744.[18] J.H. Ferziger, Computational Methods for Fluid Dynamics: J.H. Ferziger, M. Períc, Springer-Verlag, Berlin, 1996.[19] T. Maríc, J. Höpken, K.G. Mooney, The OpenFOAM Technology Primer, Zenodo, 2021, http://dx.doi.org/10.5281/zenodo.4630596.[20] L. Cao, Mesoscopic-scale simulation of pore evolution during laser powder bed fusion process, Comput. Mater. Sci. 179 (2020) http://dx.doi.org/10.1016/j.commatsci.2020.109686.[21] M. Afrasiabi, C. Lüthi, M. Bambach, K. Wegener, Multi-resolution SPH simulation of a laser powder bed fusion additive manufacturing process, Appl. Sci.(Switzerland) 11 (2021) http://dx.doi.org/10.3390/app11072962.[22] A. Queva, G. Guillemot, C. Moriconi, C. Metton, M. Bellet, Numerical study of the impact of vaporisation on melt pool dynamics in laser powder bedfusion - application to in718 and Ti–6Al–4V, Addit. Manuf. 35 (2020) http://dx.doi.org/10.1016/j.addma.2020.101249.[23] J.-P. Fürstenau, H. Wessels, C. Weißenfels, P. Wriggers, Generating virtual process maps of SLM using powder-scale SPH simulations, Comput. Part. Mech.7 (2020) 655–677, http://dx.doi.org/10.1007/s40571-019-00296-3.[24] J. Katagiri, M. Kusano, S. Minamoto, H. Kitano, K. Daimaru, M. Tsujii, M. Watanabe, Melt pool shape evaluation by single-track experiments and finite-element thermal analysis: Balling and lack of fusion riteria for generating process window of inconel738lc, Materials 16 (2023) http://dx.doi.org/10.3390/ma16041729.[25] L. Wang, Y. Zhang, W. Yan, Evaporation model for keyhole dynamics during additive manufacturing of metal, Phys. Rev. A 14 (2020) http://dx.doi.org/10.1103/PhysRevApplied.14.064039.[26] R. Murayama, E. Ohmura, I. Miyamoto, H. Hara, N. Inoue, Thermohydrodynamics analysis on mechanism of bump formation in laser texturing, Q. J. Jpn.Weld. Soc. 18 (3) (2000) 373–380, http://dx.doi.org/10.2207/qjjws.18.373.[27] H. Chen, W. Yan, Spattering and denudation in laser powder bed fusion process: Multiphase flow modelling, Acta Mater. 196 (2020) 154–167,http://dx.doi.org/10.1016/j.actamat.2020.06.033.[28] B. Mahmoodi, S. Hosseini, G. Ahmadi, CFD–DEM simulation of a pseudo two-dimensional spouted bed comprising coarse particles, Particuology 43 (2019)171–180, http://dx.doi.org/10.1016/j.partic.2017.12.014.[29] I. Mema, J.T. Padding, Fluidization of elongated particles: effect of multi-particle correlations for drag, lift, and torque in CFD-DEM, AIChE J. 67 (5)(2021) e17157, http://dx.doi.org/10.1002/aic.17157.[30] T. Yu, J. Zhao, Semi-coupled resolved CFD–DEM simulation of powder based selective laser melting for additive manufacturing, Comput. Methods Appl.Mech. Engrg. 377 (2021) http://dx.doi.org/10.1016/j.cma.2021.113707.[31] T. Yu, J. Zhao, Quantitative simulation of selective laser melting of metals enabled by new high-fidelity multiphase, multiphysics computational tool,Comput. Methods Appl. Mech. Engrg. 399 (2022) http://dx.doi.org/10.1016/j.cma.2022.115422.[32] S.L. Fuchs, P.M. Praegla, C.J. Cyron, W.A. Wall, C. Meier, A versatile SPH modeling framework for coupled microfluid-powder dynamics in additivemanufacturing: binder jetting, material jetting, directed energy deposition and powder bed fusion, Eng. Comput. 38 (2022) 4853–4877, http://dx.doi.org/10.1007/s00366-022-01724-4.[33] P.C. Carman, Fluid flow through granular beds, Trans. Inst. Chem. Eng. 15 (1937) 150–166.[34] S. Whitaker, The Forchheimer equation: A theoretical development, Transp. Porous Media 25 (1996) 27–61, http://dx.doi.org/10.1007/BF00141261.[35] S. Ergun, Fluid flow through packed columns, Chem. Eng. Prog. 48 (1952) 89–94, http://dx.doi.org/10.1029/JB088iS01p0B353.[36] M. Yang, L. Wang, W. Yan, Phase-field modeling of grain evolutions in additive manufacturing from nucleation, growth, to coarsening, npj Comput. Mater.7 (2021) http://dx.doi.org/10.1038/s41524-021-00524-6.[37] Th. Camus, D. Maisonnette, O. Baulin, O. Senninger, G. Guillemot, Ch.-A. Gandin, Three-dimensional modeling of solidification grain structures generatedby laser powder bed fusion, Materialia 30 (2023) 101804, http://dx.doi.org/10.1016/j.mtla.2023.101804.14http://dx.doi.org/10.1016/j.powtec.2021.11.056http://dx.doi.org/10.1016/j.jmapro.2023.03.040http://dx.doi.org/10.1016/j.jmapro.2023.03.040http://dx.doi.org/10.1016/j.jmapro.2023.03.040http://refhub.elsevier.com/S2214-157X(23)00783-9/sb15http://refhub.elsevier.com/S2214-157X(23)00783-9/sb16http://dx.doi.org/10.1063/1.168744http://refhub.elsevier.com/S2214-157X(23)00783-9/sb18http://dx.doi.org/10.5281/zenodo.4630596http://dx.doi.org/10.1016/j.commatsci.2020.109686http://dx.doi.org/10.1016/j.commatsci.2020.109686http://dx.doi.org/10.1016/j.commatsci.2020.109686http://dx.doi.org/10.3390/app11072962http://dx.doi.org/10.1016/j.addma.2020.101249http://dx.doi.org/10.1007/s40571-019-00296-3http://dx.doi.org/10.3390/ma16041729http://dx.doi.org/10.3390/ma16041729http://dx.doi.org/10.3390/ma16041729http://dx.doi.org/10.1103/PhysRevApplied.14.064039http://dx.doi.org/10.1103/PhysRevApplied.14.064039http://dx.doi.org/10.1103/PhysRevApplied.14.064039http://dx.doi.org/10.2207/qjjws.18.373http://dx.doi.org/10.1016/j.actamat.2020.06.033http://dx.doi.org/10.1016/j.partic.2017.12.014http://dx.doi.org/10.1002/aic.17157http://dx.doi.org/10.1016/j.cma.2021.113707http://dx.doi.org/10.1016/j.cma.2022.115422http://dx.doi.org/10.1007/s00366-022-01724-4http://dx.doi.org/10.1007/s00366-022-01724-4http://dx.doi.org/10.1007/s00366-022-01724-4http://refhub.elsevier.com/S2214-157X(23)00783-9/sb33http://dx.doi.org/10.1007/BF00141261http://dx.doi.org/10.1029/JB088iS01p0B353http://dx.doi.org/10.1038/s41524-021-00524-6http://dx.doi.org/10.1016/j.mtla.2023.101804 Influence of recoil pressure, mushy zone flow resistance and reflectivity on melt pool shape in laser powder bed fusion simulation Introduction Materials and Methods Discrete element modelling of powder bed formation Governing equations and physical models Physical properties and calculation parameters used Simulation Single-track experiment Results and Discussions Case study 1: Influence of recoil pressure Case study 2: Influence of mushy zone flow resistance Case study 3: Influence of laser reflectivity Conclusions CRediT authorship contribution statement Declaration of competing interest Data availability Acknowledgements References