LaserbeamFoam: Laser Ray-Tracing and Thermally Induced State Transition Simulation Toolkit

The application of high energy density photonic sources to the surface of metallic substrates causes localised topological evolution as the interface deforms due to hydrodynamic forces through fusion and vapourisation state transitions. Understanding how this laser energy is deposited, which may involve multiple reflection events, coupled with a thermal-fluid-dynamics framework capable of describing the heat and mass transfer in the system, permits accurate predictions of many important processes, including Laser Powder Bed Fusion, selective laser melting and laser welding among many others. In this work


Motivation and significance
High energy density electrical arc, electron beam and laser energy sources are used in various advanced manufacturing and joining processes [1].Developing high-fidelity mathematical frameworks that describe how these sources induce state transitions and substrate evolution is crucial to understanding how processing routes determine the final thermo-mechanical, topological, and microstructural state of the metallic substrate.Recently, for electromagnetically driven arc welding processes (for which there exist large gradients in the electromagnetic properties at the interface between the substrate and shielding gas), such a high fidelity multi-component magneto-thermal-hydrodynamics description has been proposed and implemented numerically [2].This permits the arc-welding heat source and momentum contributions to be explicitly captured in a fundamental manner.
For electron beam processes, a surface Gaussian heat source description, often with a Monte-Carlo model for the electron beam energy deposition, is accepted as the highest fidelity method to represent the incident electron beam at the continuum scale [3][4][5].While for laser-based processing, the highest fidelity method to numerically describe the laser source is to discretise the heat source into a finite number of rays and use a ray-tracing algorithm [6][7][8][9][10].
High energy density laser heat sources are used in various manufacturing applications, including -but not limited to -laserpowder bed fusion, selective laser melting, laser-directed energy deposition, laser welding, and laser drilling processes [1,6,11].The high energy density laser source causes localised melting, and with sufficient power density, vapourisation of the metallic substrate [12].The vapourisation of the substrate can lead to the formation of a thermo-capillary, or 'keyhole', which is the primary mechanism for the extensive penetration of welds produced using laser and electron beam sources.In the fluid state, the dynamics are dominated by the interfacial forces between the molten liquid and the vapour or gaseous phase [11].These interfacial forces are the surface tension and its dependence on temperature, which generates strong surface shear flows, and the recoil momentum induced as portions of the metallic substrate vapourise [6,[11][12][13].Other hydrodynamic forces, such as buoyancy, are present; however, their contribution to the momentum field is much smaller.
When a beam of photons is incident on an opaque substrate, a fraction of these photons are absorbed, and the rest reflected.The absorptivity describes the fraction of the beam energy absorbed by the material [14,15].The absorptivity for a given wavelength of incident radiation is a complex function of the incidence angle between the beam and the substrate and fundamental material parameters of the substrate [14].This is further complicated as the topology of the substrate will change through state transitions and as a result of the associated hydrodynamic forces.As the topology of the substrate evolves, so does the local incidence angle between the beam and substrate.Therefore the (often multiple) reflections, and associated energy deposition, that the incident laser beam experiences is a highly complex problem.This is especially important in keyhole welding scenarios where hundreds of reflections inside the thermo-capillary may occur before the packet of laser energy is either fully absorbed or leaves the keyhole structure.These multiple reflections and constantly evolving topology make the accurate representation of the laser heat source and, therefore, the prediction of the heat and mass transfer in processes driven by laser-substrate interactions challenging.Typically surface Gaussian distributions are used to represent the laser heat source.Recently these Gaussian distributions have been extended to account for the angle of incidence and the effect of absorptivity [14].
In this work, we use a ray-tracing approach and discretise the laser beam into a finite number of 'rays'.These rays are then tracked through all their reflections inside the computational domain, and the heat deposited into the substrate is computed discretely.The absorptivity of each ray at each interaction event is calculated using the Fresnel equations.The thermal-fluiddynamics framework, and absorptivity model, will now be briefly described.A software description is provided, and a series of representative examples in two and three dimensions is presented.Finally, the impact and conclusions of the work are presented.

Model description
The thermal-fluid dynamics portion of laserbeamFoam concerns the solution of equations that describe the conservation of momentum, mass and energy (namely the Navier-Stokes equation, a volume of fluid equation and an energy transport equation) in a coupled fashion.This work simplifies the vapourisation state transition to utilise a phenomenological recoil pressure term.These equations and their solution algorithm are fully described elsewhere [12]; however, an overview of the equations will now be presented for completeness.Eq. ( 1) describes the conservation of linear momentum in the system [16,17].
where U is the velocity, ρ is the mass density, P is the fluid pressure, and τ is the viscous stress tensor.Various additional source terms are present in Eq. (1) to capture important physics.S m is the Carmen-Kozeny sink term that reduces the velocity to zero in the cells filled with solid.F g describes buoyancy forces through a Boussinesq approximation.F s describe surface forces (surface tension and its dependence on temperature and recoil forces due to vapourisation) and is given as: where σ is the surface tension, κ is the curvature, P v is the recoil pressure, n is the interface normal and dσ dT is the surface tension variance with temperature.It is assumed that all phases (solid, liquid and gas) are incompressible, so a divergence-free stipulation on the velocity field is used for the system's closure, as shown in Eq. (3).
The volume of fluid approach is adopted to track the metallic interface using a volume fraction-based indicator function, α, which represents the volume fraction of the metallic phase: Additionally in Eq. ( 4) a compressive interface sharpening term is included to counteract the effect of numerical diffusion of the interface and is discussed in detail elsewhere [18].Finally, Eq. ( 5) describes the conservation of energy in the computational domain: where c p and k are the specific heat and thermal conductivity of the multi-phase mixture respectively.S h captures the latent heat effects of fusion and vapourisation state transitions in the substrate.q is the heat source used to describe the heat input to the metallic surface in the process.Typically, the energy profile across a laser beam is Gaussian in nature.The initial energy of a given ray in the beam, q, is given by: where r is the radius of the Gaussian distribution and Q is the total beam power.∆ is the computational cell length scale where the heat is deposited (1/∆ is equivalent to |∇α| in the thin interface limit).The heat source travels along the positive z-axis at a velocity v and is offset initially by distance b g and l g in the x-axis and z-axis, respectively.Σ is a parameter that is used to define the beam radius.If the beam radius is defined as the location at which the power falls to 1/e 2 then Σ = 2. Sometimes it is convenient to define the beam radius as the location at which the power falls to 5% of the peak value, in this case Σ = 3.The value of Σ can be specified through the Radius_Flavour variable in the laser properties dictionary.

Ray-tracing and absorptivity implementation
The fraction of the energy absorbed from a particular ray by the substrate, the absorptivity a, in a particular region is a complex function of the fundamental properties of the substrate and the angle of incidence of the laser heat source [7,8].
In this work, the Gaussian laser heat source, described mathematically in Eq. ( 6), is discretised into 'rays' at the boundary of the computational domain.The number of 'rays' that are propagated through the domain is directly related to the mesh-resolution, and the number of computational cells that are within the laserbeam radius at the boundary through which the laser-beam enters the domain.The rays are all given an initial vector along which they will travel through the computational domain.They are tracked along this vector until they encounter the metallic phase.A local mesh search is performed to track each ray through the domain; if the local search fails then the code reverts to a global (N 2 ) search.For each ray, upon reaching a computational cell containing the metallic phase, the reflection vector is computed using the interface normal, n, given by n = ∇α /|∇α|.
where V R is the reflection ray vector and V I is the incident ray vector.Additionally, the absorptivity is computed using Fresnel's equations -which will be described shortly -to determine the fraction of the discretised rays photons that will be absorbed and reflected along the new reflection vector respectively [7,14,15].This process is repeated until each ray in the beam has either deposited all of its energy or has gone through (often multiple) reflections and exited the computational domain.This ray-tracing of all the rays in the laser-beam is performed at every time-step in the simulation.
The absorptivity of a particular ray, a, is given by: where, according to Fresnel's reflection equations, R s and R p are the reflectance for parallel and perpendicularly polarised light defined as: where θ is the incidence angle between a particular discretised ray and the substrate where the ray interacts.χ and ψ are functions of the refractive index, n and the extinction coefficient, k, of the irradiated material [14].The form of χ, ψ, n and k is as follows: where e r and e i are the real and imaginary parts of the relative electric permittivity respectively [19], defined as: Here ω p , f and δ are the plasma frequency, incident radiation angular frequency and damping frequency, respectively and are given by:  N e is the mean number density of free electrons in the irradiated material, q e and M e are the charge and mass of an electron, respectively, ε 0 is the vacuum permittivity, c is the speed of light, λ is the wavelength of the incident radiation and R e the electrical resistivity of the irradiated substrate material.Fig. 1 shows how this absorptivity changes as a function of the angle of incidence for two common wavelengths of incident laser radiation; the material properties used to generate this figure was chosen to be representative of 316L stainless steel.
It can be seen from Fig. 1 that the absorptivity does not change significantly at small angles of incidence.Furthermore, the angle of incidence corresponding to the peak absorptivity for the substrate material varies depending on the wavelength of the incident radiation.Fig. 2 shows a laser beam and its reflection that has been discretised into 4220 rays, incident on an interface that has deformed through melting and vapourisation eventsrendering the interface non-planar.It can be seen from Fig. 2 that with once the substrate topology changes due to melting and vapourisation events, the incident rays generally each encounter a different interface normal vector where the ray meets the substrate, and therefore the reflection of each ray is different.Note, for post-processing, the laser beam paths for all discretised beams are written out in two scalar fields, namely Ray_Q and Ray_Number; these fields can be used to visualise the paths the beams took.However, these scalar fields are over-written in the order of the ray index (indicated in the Ray_Number field).Therefore the reflection behaviour of rays with a low index may often be over-written by the rays with a higher index.It is possible to write all ray paths out entirely; however, this would use significant memory (the beam may be discretised into tens of thousands of rays in specific applications).Therefore it was decided in the development that it was preferred to write out two scalar fields for the reflection behaviour and maximise the performance of the code.Fig. 3 shows an example of the Ray_Q and Ray_Number fields for a scenario where a Gaussian beam is incident to a metallic substrate at an oblique angle.
As previously stated, the beam paths are written sequentially to these two files to minimise memory and hard disk usage.As shown in Fig. 3, the paths of the rays with a higher index have over-written the paths of rays with a lower index.For completeness, Fig. 4 shows the path of a single ray in a Gaussian beam (coloured by the ray energy normalised with respect to the initial beam energy in red) as it undergoes multiple reflections inside a thermo-capillary before finally exiting the keyhole (blue).
As can be seen from Fig. 4, each time the discretised ray reflects at the metallic interface, a portion of its energy is deposited into the substrate, and the path taken by the ray is a complex one with many reflections before the ray finally exits the domain; the ray exiting the keyhole is coloured blue representing the lower energy state.It would not be computationally tractable to write all ray paths out in the manner shown for a single ray in Fig. 4; hence the amalgamation of all ray paths into two output files.

Software description
laserbeamFoam is free software distributed under the GNU general public license.laserbeamFoam has been developed using the OpenFOAM v6 libraries; however, a branch is available in the repository that compiles against the newer OpenFOAM-10 (and the ESI version) libraries.laserbeamFoam can be compiled the OpenFOAM wmake build command.Full installation instructions are provided in the repository.The code is intended for use in applications where the thermal-fluid dynamics of processes driven by laser-substrate interaction are of interest, such as Laser-Powder Bed Fusion or laser drilling applications.However, the software could be used for many applications beyond these areas where laser-substrate interaction and state transitions are of interest.Within the repository are a series of illustrative examples; it is intended that these examples will grow in number through community additions and contributions by the development team.Each commit to the laserbeamFoam repository is subject to a code review, mathematics review, and regression testing in the development branch before it is added to the main branch in the GitHub repository, ensuring code stability between versions.laserbeamFoam is fully parallelisable using the message passing interface (MPI) framework.
It is expected that a laserbeamFoam user will wish to apply the framework to their problem of interest.The easiest way to do this is to copy an example to a new directory, modify the domain, material, and laser description files, and then run the laserbeamFoam solver on the new case.To perform a simulation using laserbeamFoam (after the solver is compiled against the OpenFOAM libraries), the following steps must be taken: 1. Navigate to a model folder, e.g.tutorials/plate3D.$> setFields 6. Run the solver: • For parallel runs: -Decompose the problem into the required number of sub-domains (only required for parallel runs): $> decomposePar -Execute the solver in parallel with MPI (here six processes are used) $> mpirun -np 6 laserbeamFoam -parallel >log & • For serial runs: -Execute the solver in serial $> laserbeamFoam 7. Outputs will be written in the same directory.
Note that the boundary and initial conditions for the temperature, metallic volume fraction (α), velocity and pressure fields can all be edited in the initial files before any simulation run.Furthermore, the boundaries at which the laser source enters the domain are set in the initial folder by setting the Laser_boundary field to a positive number at the boundaries at which the laser enters the domain.Further details on the meshing, decomposition, and field operations can be found in the OpenFOAM documentation.The OpenFOAM framework was chosen as it is open-source, highly parallelised, widely used and validated in academia and industry.The vibrant community of researchers active on web forums like www.cfd-online.com is an additional incentive to use the OpenFOAM framework.
laserbeamFoam is fully parallelised, and OpenFOAM offers several options to decompose the simulation domain into subdomains such as simple, hierarchical, metis and manual.There is no limit on the decomposition chosen for the parallelisation; the ray-tracing algorithm has been developed to work with all decomposition schemes.

Illustrative examples
In the following examples, all material properties represent 316L stainless steel substrate under an Argon gas atmosphere.

2D plate
In this example, a two-dimensional domain is used to investigate the keyhole formation and evolution as a 2 kW laser with wavelength 1.064 × 10 −6 m and beam radius of 5.0 × 10 −4 m is incident on a 3.0 × 10 −3 m thick substrate.In the first instance, the discretised Gaussian beam is incident normal to the substrate interface; a subsequent simulation shows the evolution behaviour with the same beam incident at 45 • to the interface.Fig. 5 shows the simulation domain at three time-steps as the beam forms a thermo-capillary and penetrates through the domain.
The energy from the laser heat source is deposited in a complex manner inside the thermo-capillary due to the many reflections within the keyhole structure.Fig. 5 shows that the laser heat source penetrates the domain in 0.03 s.In the second instance, the angle of incidence of the beam is changed to an oblique angle.Fig. 6 shows the simulation domain at three time-steps as the beam forms a thermo-capillary and penetrates through the domain.
Fig. 6 shows that the beam penetrates the domain in the oblique scenario in 0.047 s.As can be seen through comparison between Figs. 5 and 6, when the initial beam incidence angle is changed from 0 • to 45 • the overall penetration time increases by approximately 50%.The simulations help to elucidate the complex reasons for this.In the oblique scenario, the Gaussian heat source is spread over a larger area, so the energy density, in this case, is reduced relative to the 0 • scenario.Furthermore, in the reference frame of the beam, the substrate is thicker in the oblique case.The fluid dynamics of the molten substrate also play an important role in determining the penetration rate (the penetration times are not linear functions of substrate thickness in the beam direction); molten material on the side-walls of the thermo-capillary is more likely to flow back into the keyhole in the oblique case -decreasing the penetration rate.Animations of these scenarios are included in this work.

3D plate
In this scenario, the solver is applied to simulate the evolution of a three-dimensional substrate when a 5 × 10 3 W beam is directed at the plate at a 45 • angle of incidence.The substrate is 3×10 −3 m thick and is surrounded by an Ar atmosphere.The laser heat source is static and discretised into approximately 4600 rays (this figure is determined from the mesh density).Fig. 7 shows snap-shots of the evolution of the temperature and topology of the substrate.
As shown in Fig. 7, applying the laser heat source at the plate surface causes a thermo-capillary to form.The oblique angle of incidence of the beam to the plate means that the dynamics of the generated keyhole are highly irregular, as was observed in the 2D case described previously.

2D simple powder particles
In this example, the laser heat source passes over the substrate at an oblique angle in a two-dimensional domain.The domain contains seven circular metallic substrate particles with two different radii on top of a planar metallic substrate.Fig. 8 shows the domain at t = 0 s and five subsequent time-steps as the powder particles melt due to the traversing laser heat source.
As can be seen from Fig. 8, the incident discretised laser heat source deposits its energy into the substrate and initiates melting and vapourisation.There are often complex reflection paths as the individual rays may traverse between the circular particles and undergo many reflections, therefore complex absorptivity and energy deposition mechanisms.The reflection paths of the discretised rays are also shown in Fig. 8; however, as previously discussed, the higher index rays obscure the reflection characteristics of the lower index rays.An animation of this example application accompanies this manuscript that shows the porosity that becomes trapped between the particles escaping.

Laser-powder bed fusion
In this example, application laserbeamFoam simulates a laserpowder bed fusion process.Here a 500 W laser source is discretised into 300 discrete rays.The heat beam travels along the domain at 1.0 ms −1 .Fig. 9 shows the evolution within the computational domain.
As shown in Fig. 9, the effect of Marangoni flow is pronounced; surface shear flows are generated due to the differential heating along the metallic-gaseous interface.This has the effect of pushing metallic material to the right as time progresses.It is also interesting to note that it is possible to predict porosity and lack of fusion, as seen in Fig. 9.

Validation against in-situ laser welding data
In this final illustrative example, we will validate the proposed software against a well documented laser welding experiment where the keyhole morphology and penetration rate was observed using ultrahigh-speed synchotron X-ray imaging.In this case a 3D simulation was utilised, and the beam radius was defined as the distance at which the peak intensity decreased to 1/e 2 of the peak value.In the experiment a 156 W laser beam with a spot-size of 140 µm was incident on a Ti6Al4V substrate.
The full experimental set-up is described elsewhere [20].Fig. 10 shows the numerically predicted, and experimentally observed, thermo-capillary morphologies at a number of instances.
It can be seen from Fig. 10 that the proposed solver accurately predicts the thermocapillary morphology and penetration rates, particularly the change in penetration rate as more of the rays from the beam are reflected towards the base of the keyhole.The solver also predicts the change in oscillation modes of the keyhole  for the volumetric dilation and vapour plume generation.For a more accurate prediction of these systems, that are dominated by vapourisation and condensation events, the ray-tracing algorithm will be incorporated into a more robust mathematical treatment of the vapourisation and condensation state transitions (where the volumetric dilation due to these events is explicitly included), as was presented elsewhere [11].

Impact and conclusions
The presented open-source thermal fluid dynamics solver, incorporating a complex ray-tracing description of the heat source deposition through multiple reflections, enables a wide range of researchers to explore complex laser-substrate interactions.The presented solver could be utilised in a wide range of applications in the manufacturing community, including (but not limited to): creating digital twins of the laser-powder bed fusion process through higher fidelity representations of the topological evolution and energy deposition through complex reflections, investigating the penetration rates of laser sources through metallic substrates in laser welding and drilling scenarios, predicting the porosity propensity in laser welding applications.The presented solver will be useful if the fundamental physics of laser beams interacting with alloy substrates is of interest.laserbeamFoam is part of the OpenFOAM landscape, and therefore it is expected that the solver will evolve further with the inclusion of additional physics and community suggestions; for example, in the future, the electrical resistivity, as well as other material parameters could readily be made functions of temperature.The authors are open to all suggestions from the community.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Absorptivity as a function of incidence angle, θ, for two different wavelengths of incident electromagnetic radiation.

Fig. 2 .
Fig. 2.An incident beam travelling from the upper left in the image meets a deformed interface which scatters the incident radiation in a complex manner.The discretised heat source rays are coloured by their respective energy.Note that only 10% of the rays are shown for image clarity.

Fig. 3 .
Fig. 3. Example output for the Ray_Q and Ray_Number fields from the laser-beamFoam solver; for a Gaussian beam, discretised into 96 rays, incident on a metallic substrate (shown in grey) at an oblique angle.The beam travels from the top-right in these images, interacting with he substrate and depositing a fraction of each 'rays' energy, and then reflecting and exiting the domain at the top left.

Fig. 4 .
Fig. 4. Schematic Diagram showing the route of a single discretised ray in the Gaussian beam showing an representative path inside the thermo-capillary structure in a laser welding scenario.The magenta line delimits the solid/liquid states, and the green contour delimits the metallic/non-metallic phases.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

2 . 5 .
Remove any previous simulation files $> rm -r 0* processor * 3. Copy boundary and initial condition into the initial time folder: Set the initial conditions:

Fig. 5 .
Fig. 5. Evolution of the thermal field and topology of the metallic substrate as a laser heat source is applied normal to the substrate interface.Temperature is shown in Kelvin.The multiple reflections of the discretised laser source can be seen, coloured by the energy of the individual beams, showing the complex laser-substrate interaction that produces the thermo-capillary.The magenta line in the sub-figures shows the ϵ 1 = 0.5 contour (delimiting the solid and fluid states), while the white line shows the α = 0.5 contour (delimiting the metallic and non-metallic phases).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 6 .
Fig. 6.With the initial angle of incidence of the beam changed to 45 • the keyhole structure, and penetration rate are markedly different to the case with 0 • incidence angle.The magenta line in the sub-figures shows the ϵ 1 = 0.5 contour (delimiting the solid and fluid states), while the white line shows the α = 0.5 contour (delimiting the metallic and non-metallic phases).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 7 .
Fig. 7. Time evolution of a three-dimensional thermo-capillary generated through the application of a 5 × 10 3 W laser with a 45 • angle of incidence.

Fig. 8 .
Fig. 8. Example application of laserbeamFoam where a Gaussian laser source passes over circular regions of varying diameter representative of powder particles in the powder bed fusion process at various time-steps.The laser source is incident on the substrate at an oblique angle.The magenta line in the sub-figures shows the ϵ 1 = 0.5 contour (delimiting the solid and fluid states), while the white line shows the α = 0.5 contour (delimiting the metallic and non-metallic phases).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 9 .
Fig. 9. Evolution of a series of 316L stainless steel powder particles as a laser heat source traverses the domain and causes state transitions and topological changes.The magenta line in the sub-figures delimits the solid and fluid regions, and the white lines delimit the metallic and shielding gas phases.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)that were observed experimentally.However, the framework appears to slightly over-predict the width of the thermo-capillary; This over-prediction of the width is because the treatment of the vapourisation state transition is relatively simple in the opensource version of the solver and does not explicitly account

Fig. 10 .
Fig. 10.Example application of laserbeamFoam solver used to simulate an in-situ laser generated thermo-capillary.