An extended laser beam heating model for a numerical platform to simulate multi-material selective laser melting

A laser beam heating model (LBHM) is an important part of a platform for numerical modelling of a multi-material selective laser melting process. The LBHM is utilised as a ray-tracing algorithm that is widely applied for rendering in different applications, mainly for visualisation and very recently for laser heating models in selective laser melting. The model presented in this paper was further extended to transparent and translucent materials, including materials where transparency is dependent on the material temperature. In addition to reflection and surface absorption, commonly considered in such models, phenomena such as refraction, scattering and volume absorption were also implemented. Considering associated energy transfer, the model represents a laser beam as a stream of moving particles, i.e. photons of the same energy. When the photons meet a boundary between materials, they are reflected, absorbed or transmitted according to geometric and thermal interfacial characteristics. This paper describes the LBHM in detail, its verification and validation, and also presents several simulation examples of the entire selective laser melting process with implemented LBHM.

Probability of reflection or transmission p d and p g Probabilities of the diffuse and glossy reflection or transmission (p d + p g = 1) p 1 and p 2 Uniformly distributed random numbers in the range 0 to 1 R Direction of specular reflection r Radius of the sphere V Outgoing direction x Distance from the vertical axis of the sphere, coordinate x' Relative distance from the vertical axis.

Special characters α Attenuation or absorption coefficient ϑ
Zenith angle, the angle between the specular R and outgoing V directions ρ Reflection coefficient τ Transmission coefficient φ Azimuth angle Ω Solid angle Subscripts a Absorbed radiance BRDF Bidirectional reflectance distribution function

Introduction
Many material properties such as heat conductivity, heat capacity, density and surface tension are temperature dependent, making the temperature and temperature gradients the most crucial process parameters in selective laser melting (SLM). These temperature-dependent properties influence and determine the thermomechanical behaviour, hydrodynamic effects and finally the quality of the final product. Due to such a large impact, mathematical modelling of thermal interactions during SLM is significantly important [1]. It is difficult to overestimate the importance of mathematical modelling in understanding and prediction of the complex relationship between different physical phenomena associated with this process. The physical phenomena involved and their influence on the final microstructure and properties of the metallic components manufactured by SLM can be found in [2]. Offering a unique geometrical flexibility for designers, the SLM process parameters can be altered to design an appropriate local microstructure and custom-made parts to individual requirements. Process planning by using conventional experience-driven approaches is sometimes very expensive, prompting development of computational simulation techniques to save time and resources.
A variety of numerical methods and techniques, such as finite volumes (FV), Lattice-Boltzmann (LBM) and arbitrary Lagrangian-Eulerian (ALE), have been applied to tackle different challenges in modelling of SLM and similar additive manufacturing processes at the powder scale associated mainly with free surface evolution and surface tension [3][4][5].
Powder deposition with subsequent laser treatment has also been investigated using the discrete element method (DEM) [6][7][8]. The DEM method allowed for taking into account the statistical nature of the powder bed. However, a physical representation of the melt flow was not accounted for. Hybrid DEM-meshfree continuum-based numerical models were developed to tackle the mentioned drawback, proving the feasibility of such approaches for modelling the fusion of metal powder particles in the additive manufacturing processes of this type, highlighting at the same time the importance of the appropriate heat source modelling. Lasers can generate power continuously and in pulses [9]. This study is focused on modelling the continuous laser beam as a concentrated heat source and thermal interactions between such a concentrated heat source and the generated powder bed during the SLM process. In many modelling cases, mainly due to the simplicity of their implementation, volumetric heat sources are considered based on empirical observations [10]. The significant problem in SLM modelling is related to energy absorption, which is highly dependent on the overlap of the volumetric intensity distribution with the irradiated geometry [11], presenting significant obstacles in powder scale SLM simulations. One of the solutions to improve accuracy is to describe the laser beam by Maxwell equations representing its electro-magnetic field components, as has been shown elsewhere [12]. Although this formulation allows for achievement of a high level of accuracy, the approach cannot be considered acceptable for engineering practice due to the high computational cost. For a more accurate numerical consideration of absorption and reflection in SLM, ray tracing is the most appropriate method where the beam is discretised into portions of energy called rays. Assuming that diffraction is negligible, ray tracing is a purely geometric approach widely used in optics [13][14][15] and also for investigating incident radiation scattering in particulate domains [16][17][18]. In metal processing simulations, a ray-tracing algorithm has been applied in DEM to model the interactions between the laser beam and the irradiated powdered material [19]. Another example is the model developed for the investigation of laser keyhole welding process based on the finite difference method, where the evolution of the liquidvapour interface is tracked with a level-set method providing the surface normals by using a ray-tracing algorithm [20]. The ray-tracing method has been applied to determine the absorptivity of metal in laser welding, which was then used in numerical modelling based on the smoothed particle hydrodynamics (SPH) method [21,22]. In many SLM processes, the wavelength of the incident laser radiation is more than an order of magnitude lower than the diameter of powder particles and relevant features of the surface that make a ray-tracing approach appropriate for physically based numerical modelling.
Laser heating models are a very important aspect in the numerical simulation of SLM. One common and efficient technique is the estimation of the absorption coefficient following the ray-tracing model (RTM). Conceptually, in the RTM, the laser beam is divided into rays with a certain size, direction and amount of power, incident on the surface of the powder bed [23]. Yang et al. [24] proposed a 3D optical model, using the ray-tracing algorithm to simulate the energy interaction of the powder bed particles and the laser beam to investigate the laser energy absorption behaviour of the SLM-processed AlSi12 powder. Le et al. [25] presented a model considering the interactions between the electron beam and Ti6Al4V based on Monte Carlo ray-tracing simulations and the total absorption of Ti6Al4V material on microscale. At the mesoscopic level, the absorptivity profile along the depth of the powder bed was calculated using the absorptivity of Ti6Al4V for various incidence angles and the 3-D powder bed model. In both articles, a higher absorptivity of the powder is observed compared to the bulk material [24,25]. Liu et al. [26] presented a ray-tracing heat source model used in mesoscale finite volume method (FVM) simulations for the selective laser melting (SLM) process. 3-D mesoscale simulations were performed for the 316L stainless steel.
Zhou et al. [27] analysed the discrepancies between the Gaussian surface heat source model and the ray-tracing heat source model for numerical simulation of selective laser melting. They show that structures with more sheltered surfaces (by adjacent structures) or higher degrees of self-shielding will lead to larger discrepancies between the simulation results with these two heat source models. For the single metal sheet, there are few discrepancies. For a powder bed with complex particle distribution, the Gaussian surface heat source model adds more heat to the powder bed compared to the ray-tracing heat source model, while the discrepancies decrease gradually as the metal powder melts and the fluid fills the gaps. It should be noted that the publications [19][20][21][22][23][24][25][26][27] consider the RTM for application to opaque materials with reflection, but not to transparent and translucent materials.
Liu et al. [28] predicted the porosity in SLM by a combination of experimental, numerical and analytical calculations. A CFD modelling was performed to simulate the single-track scanning to reveal the physical mechanisms of the defect formation in the SLM process. The single-track SLM tests were simulated on a FLOW3D finite volume method (FVM) platform by solving the Navier-Stokes and heat conduction equations.
Sagar et al. [29] presented a simulation method on the simulation platform based on DEM, CFD and structural mechanics. Packing density, powder layer thickness, melt pool layer thickness, layer displacement and the solidified layer thickness were measured as responses from the simulation study. The combination of DEM for powder spread and CFD for melting/solidification models was presented by Li et al. [30]. The sensitivity of modelling terms such as melting and solidification mass transfer coefficients and the mushy zone constant is tested, and their effects are quantified in terms of melt pool dynamics, solidified surface morphology and pores formation.
Fotovvati et al. [31] presented a study of the influence of SLM parameters on the surface characteristics of 'singlelayer' Ti6Al4V powder using experimental, computational and data-driven approaches. The computational approach uses the DEM for the powder bed distribution and the CFD for simulating the laser-powder interaction. Fotovvati and Chou [32] presented a multilayer model of the powder bed fusion process. The approach uses DEM and CFD in a sequential manner to about 10 layers.
Recently, multi-material AM processes were also modelled. Shinjo and Panwisawas [33] using the coupled levelset and volume-of-fluid method for multiphase thermal fluid flow simulation of very complicated AM processes simulated and analysed the fluid flow dynamics and mass loss characteristics due to metal vaporisation. They obtained, among others, the results connected with the formation of the melt pool, keyhole and flow for different elements: the mass loss rate as a result of vaporisation. A multi-material model for the lattice-Boltzmann-based simulation of the mixing of two elements to a binary alloy via PBF is presented by Küng et al. [34]. Tang et al. [35] simulated the multi-material laser PBF process using the volume of fluid CFD model. They modelled two-, three-and four-component powders with mixed titanium Ti, niobium Nb, vanadium V and arsenic Ar particles. In real SLM processes, materials can be mixed uniformly before the process, especially when they are metals and can form alloys. When materials such as metal and glass (ceramics) are used, where materials have very different properties and should not form alloys, and the structure of the product should be accurately controlled because materials play different roles in the product, the mixture before the process cannot be performed [36].
Recent publications [28][29][30][31][32][33][34][35][36] presented mainly starting from 2022 testify to a new stage in SLM modelling and indicate an exit to a new level. Our team started moving in the same direction in 2017 [36,37] and we believe that we are now at the same stage and at the same level [38,39].
The objective of this paper is a presentation of the RTM-based light propagation model implemented into the 1 3 platform for modelling the SLM process. RTM is developed to investigate the thermal interactions between a laser beam and a powder bed during SLM, which are characterised by the melting and solidification of the deposited particulate material. Considering associated energy transfer, the model represents the laser beam as a stream of moving particles, i.e. photons of the same energy interacting with a deposited powder bed in the liquid or solid states. The other behaviour of the beam, such as reflection, absorption or transmission at the interface, is predicted by assuming the appropriate material properties, including geometric and thermal interfacial characteristics. The proposed model is exemplary and being imbedded into the multiphysics model of an entire powder bed fusion process based on homogeneous cellular automata (CA) and lattice-Boltzmann (LBM) methods represents a viable numerical tool for prediction of thermal interactions during SLM simulations. The outline of this work is as follows. Section 2 presents an overview of SLM modelling based on CA and LBM methods. The model of light propagation is discussed in detail in Section 3. Several verification methods are discussed in Section 4 followed by model validation (Section 5) based on comparison of visual and numerical picture analysis, theoretical explanation and simulation results. Section 6 illustrates an example of the modelling of the entire SLM process after implementation of the discussed model of light propagation into the holistic model of the considered powder bed fusion process. The article concludes with an outlook and some recommendations for further development (Section 7).

A role of the laser beam heating model in the modelling of selective laser melting
The paper discusses the development of the laser beam model and its thermal action on powder in the context of the multiscale numerical modelling of the SLM processes characterised by melting and solidification. Several physical events are considered in the holistic model of the SLM process, including powder bed deposition, laser energy absorption and heating of the powder bed by the moving laser beam leading to powder melting, fluid flow in the melted pool, flow through partly or not melted materials, and solidification. The principles of the holistic model of the SLM process are discussed in our previous publication [36,37,40]. Figure 1 represents the proposed diagram, where the process is divided into different stages according to the associated physical phenomena, which are related to the corresponding mathematical models that constitute the holistic process model. Among other elements, the presented diagram includes items marked by green colour, which are the object of this paper, namely the laser beam treatment, which is linked to the associated physical process and the relevant mathematical model. This study focuses on the development of the highlighted laser beam heating model that presents relevant simulation examples within the framework of the complete numerical platform.
The modelling approach is based mainly on two homogeneous methods, CA and LBM. Components of the model operate in the common domain, allowing them to be linked into a more complex holistic numerical model with the possibility of complete full-scale calculations avoiding complicated interfaces. According to these assumptions and conception, the relevant (appropriate) model of laser beam heating should be developed.
The simplified model of heat transfer, initially used in the holistic model for illustrating purposes, has been discussed earlier [23]. Particle or photon tracking methods, known as well as ray-tracing, are widely applied for visualisation and computer graphics. A similar algorithm is used in the presented model of laser beam heating. It is assumed that the laser radiation transmits light energy in gas as a straight beam without scattering. When the surface of the liquid or solid phase is struck, the energy is divided into three parts: reflected, absorbed and transmitted. If the material is opaque, the beam is not transmitted, and the radiation is reflected and absorbed. The light energy that penetrates the material is subject to scattering and absorption. Both transmission and absorption depend on a variety of chemical and physical characteristics of the powder bed material and also on the wavelength of the incident radiation.

Model description
In the model, the laser beam is considered as a stream of moving particles, photons of the same energy. It is assumed that in air (gas) the particles move straight without any scattering or aberration. When a particle meets a material in the liquid or solid state, its further movement depends on the material and surface properties. The particle can be reflected, absorbed or transmitted at the interface (boundary). Reflection, absorption or transmission is randomly chosen with probabilities depending on the appropriate physical characteristics of the material.
The bidirectional scattering distribution function (BSDF) is used to define the direction of reflection and transmission of the particle. In this study, BSDF is considered two independent functions, that is, bidirectional reflectance distribution function (BRDF) and bidirectional transmittance distribution function (BTDF) for the reflected and transmitted streams, respectively. BRDF and BTDF are functions of four real variables. The variables define the incoming and outgoing directions of light. Each direction is described by two angles, for example, the zenith ϑ and azimuth φ angles. The zenith angle ϑ is the angle between the stream direction and the normal to the surface, while the azimuth angle φ defines an angle on the surface (or tangential to the surface) between the projections of the directions of the incoming and outgoing streams. The BRDF and BTDF define how light is reflected or transmitted. The functions return the ratio of the reflected light exiting along the outgoing direction to the incident light.
The applied algorithm is similar to the photon mapping algorithm developed by Jensen as an efficient alternative to pure Monte Carlo ray-tracing techniques [41,42]. In this study, the algorithm used in modern visualisation methods has been adapted to deal with the thermal interactions between a continuous laser beam and powder material. The visualisation methods are limited to one final point such as eyes observing the scene; therefore, the method defines both directions of the bidirectional functions explicitly. In terms of thermal interactions, the outgoing direction should be chosen randomly according to the BRDF or BTDF.
The model of light propagation was developed using object-orientated programming in C++ and OpenGL version 3.3 for rendering. Additional libraries used were GLFW (Graphics Library Framework) for window management and OpenGL contexts, and GLAD for loading OpenGL functions at runtime. The Microsoft Visual Studio 2022 community version was used as the integrated development environment (IDE).

The model algorithm
In the model of light propagation, the laser beam is represented as a photon beam with given dimensions of beam and energy of the photons. The number and location of the photons depend on the parameters of the laser beam. Then, every photon is considered separately according to the algorithm presented in Fig. 2.
The model is probabilistic. The types of photon interaction with the material, the direction of the reflection, refraction, scattering and the length of the photon path are defined with the use of random numbers.
The algorithm begins with the definition of photon location and its moving direction. The photon movement depends on the material through which it is moved. Initially, it is the air (gas) that is considered transparent, i.e. the photon moves straightly without scattering, absorption or aberration. If the material is opaque, for example, metal, the energy of the photon is absorbed here. In the case of translucent material, for example bioactive glass, the photon can be scattered or absorbed; then, the lengths of scattering l s and absorption l a paths are defined based on material optic properties with the use of probabilistic methods.
Then, the photon movement is considered. There are four cases of the end of this movement considered in the model. They are the space boundary, the material boundary, scattering and absorption. The choice of the case depends on the minimal value of the lengths: scattering length l s , absorption length l a , length to material boundary 1 3 l m and length to space boundary l b . Any action is applied in the case of a space boundary; the photon simply flies out of space, and its energy is lost. The photon energy is transferred to the other model as a heat source in the case of absorption. In the case of scattering, a new length l s of scattering path is defined, and the algorithm returns to the photon movement.
The last case of the material boundary is more complicated and is considered below. Generally, incident radiance (power) can be divided into three parts: If the material is opaque, the transmittance coefficient is equal to zero, τ = 0. From the definition of the mentioned coefficients, their sum is equal to the unit: The probabilities of photon reflection p r , transmission p t and absorption p a on the surface can be set equal to the appropriate coefficient: The types of photon interaction with the material boundary can be defined by using the random number with corresponding probabilities according to (3). Correspondingly, to define photon reflection, transmission or absorption, it is enough to draw a random number in the range 0 to 1 and determine in which range this number lies (e.g. 0 ÷ p r , p r ÷ p t + p a , p r + p a ÷ 1). Reflection and transmission are described by the BRDF or BTDF, which (1) L i = L r + L t + L a = ( + + )L i .
can be represented as a sum of the diffuse f d and glossy f g reflection (transmission) functions as follows: It can be noted that the integrals of the BRDF, BTDF, as well as diffuse f d and glossy f g reflection (transmission) functions, along with all possible outgoing directions equal the unit. Then, the sum of the probabilities of the choice of the surface type that can be applied for further calculations equals the unit (p d + p g = 1). To choose the type of surface, it is enough to draw a random number p in the range 0 to 1 and determine in which range this number lies (e.g. p < p d or p > p d ). After choosing the type of reflection or transmission (diffuse or glossy), the inverse BRDF and BTDF should be used to obtain the direction based on the probability. Both functions, that is, diffuse f d and glossy f g , as well as their inverse functions using probabilities, are described below.
The reflected or transmitted photon is traced further in a new direction until it reaches another surface or the boundary of the modelled space, or it is scattered or absorbed, as described above. A location of the absorbed photon is discretised to assign its energy as a heat source to an appropriate node in the mesh of the LBM model of heat transfer. It is the end of the algorithm for the considered photon.

Model equations
The surface models with appropriate distribution functions are presented in this section. Absorption and scattering models for translucent material are considered as well. Fig. 2 Block scheme of the model algorithm There are three emission laws with correspondent distribution functions described here. The emission law with uniform distribution is used in the model as the point or volumetric light source, and it can be easily verified visually. However, it is not used in interactions of the laser (photon) beam with the material surface. Lambert's emission law is applied for a matte surface, and the Gauss distribution function is applied for a glossy surface.
The schematic representation of the different surface models with a schematic representation of the distribution functions used in this study is shown in Fig. 3. The surface model is presented as the Lambertian reflectance model for the matte surface (Fig. 3c, f) and the glossy model (Fig. 3b, e). The main reflection or refraction directions on the glossy surface are represented by dotted lines in Fig. 3b and e. Similarly, the ideal reflection and refraction are shown in Fig. 3a and d.

The uniform emission distribution function
The uniform emission distribution function means that the radiance L is the same in all directions L=const. The whole emission I can be calculated as following: A solid angle Ω (5) in spherical coordinates can be defined by the azimuth angle φ and the zenith angle ϑ. Then, the differential is of the following form: It leads to the following equation: Now, it is easy to define the cumulative distribution function for both angles: To obtain the outgoing direction, the azimuth φ and zenith ϑ angles should be defined. Then, the random number p (in the range 0 to 1) is established as equal to a value of the cumulative distribution function p = D, and an appropriate angle can be received from the inverse cumulative distribution functions. For considered cases, they take the following forms: For ϑ max = π/2, it will be:

Lambert emission law
Lambert's emission law applied for the matte surface can be written in the following form: Omitting transformations similar to (7)-(11), the azimuth angle φ takes the same form (10), while the zenith angle ϑ can be calculated from the random number p according to the following equation:

Glossy emission
There are several models used for glossy surfaces, depending on the applied principles and conditions. The gloss model can be described by using a Gaussian distribution. In this work, two Gaussian distribution functions with the same (or different) dispersions are applied. The first one describes luminance depending on the possible deviation from the main direction (defined by the reflection or refraction angle) in the plane defined by the incident direction and the surface normal. The second function defines luminance depending on the deviation in the plane perpendicular to the previous one. Then, two 'zenith' angles ϑ 1 and ϑ 2 in two orthogonal planes can be obtained from the inverse cumulative Gaussian distribution. Calculation of each of these angles requires two random numbers p 1 and p 2 . Then, the appropriate zenith angle (ϑ 1 and ϑ 2 ) can be obtained as follows:

Reflection and transmission models, BRDF and DTDF
The reflection and transmission models use the BRDF and BTDF functions according to (2)-(4). These functions represent the sum of Lambert's law (13) and the glossy model with the Gaussian distribution function. The outgoing direction is calculated using equations (10), (14) and (15).

Absorption and scattering models for translucent material
The Beer-Lambert-Bouguer law is used for absorption, which can be presented in the following form: Then, the length of the photon path l defined by the random number p in the range 0 to 1 can be calculated as follows: An equation similar to (16) can be applied for scattering. Then, the length of the photon path to the collision and scattering can be calculated according to the following: The properties (l a , l s ) of the translucent materials, e.g. bioactive glass, depend highly on the temperature, which in turn changes along the path of the photon. Then, (17) and (18) should be rewritten in the following form as a condition of absorption or scattering: where ∆x -mesh step used in holistic model for spatial discretisation.

Model verification
Verification of all different predictive aspects of the holistic model is a relatively long and time-consuming process that started with the verification of the predictive abilities of powder bed formation, including effects of particle size distributions, the speed and sequence of powder delivery during deposition, indicating improvement of the prediction accuracy regarding the powder bed packing density and morphology [37]. It has been widely recognised that quantification of the packing density and its correlation to the results of the process is one of the important challenges in improving the improvement of process reliability [43]. In this work, the following methods are chosen for verification of the model of light propagation: units test, proof of correctness and integration test. The results of the light propagation model testing are presented below.

Units tests
Unit tests are the simplest ones; they only verify that the appropriate models work properly. Units represent the correspondent components of the holistic model of light propagation, which are submodels including reflection, refraction, absorption and scattering. This kind of testing demonstrates the application of the surface model with distribution functions discussed in the previous section and covers the cases illustrated in Fig. 4. Figure 4a represents reflection from the mirror surface, while the cases of reflection from glossy and matte surfaces are shown in Fig. 4b and c, respectively. For the glossy surface, the Gaussian distribution function is applied with a dispersion of 10°. On the matte surface, the photons are dissipated uniformly over the range of 0-180°. Similar results are obtained for the refraction illustrated in Fig. 4d-f. The same dispersions of 10° and 180° are applied, respectively, for the glossy surface (Fig. 4e) and the matte surface (Fig. 4f). The refraction coefficient of 1.3 has been assumed for the refraction cases presented. The modelling results illustrating different cases of multiple reflections from the circle particles are shown in Fig. 5. In these cases, it is assumed that no absorption and transmission took place, i.e. all the photons are reflected. The same reflection parameters, such as 10° dispersion for glossy surfaces (Fig. 5a and c) and 180° dispersion for the matte surface (Fig. 5b), were applied. Note that the specular reflection and transmission is much larger than the diffuse scattering for many surfaces. In these modelling cases, both the BRDF and the BTDF are considered to be relatively smooth functions, making interpolation near the specular directions easier. One potential problem can arise when the laser source is smaller than the texture features of the material surface. The resulting BRDF and BTDF models may be more specular overall than the actual material surface.
It should be mentioned that the results obtained agree with the observations of other researchers [44,45] and completely satisfy the compatibility requirements of the holistic SLM model under consideration. Figure 6 represents several results of the 3D model of scattering in translucent material with l a >> l s . Figure 4d-f illustrate refraction without consequent absorption and scattering, while Fig. 7 shows simulation results of absorption and scattering in translucent material with different relations between absorption and scattering, namely Fig. 7a represents 3l a = l s , Fig. 7b represents l a = l s and Fig. 7c represents l a = 3l s .

Correctness tests
Correctness tests are performed to prove that the submodels give the proper results for different types of surfaces. The first stage of such verification is visual, the second one is numerical. In spite of the BRDF and BTDF applied for reflection and refraction, respectively, they consider the same types of surface and use the same equations. Thus, only reflection is tested here. The screen absorbing the reflected photons is set above the flat surface. The absorbed photon changes the intensity of a colour of the appropriate pixel on the screen.
The results of the simulation presented in Fig. 8 illustrate two variants of the surface for visual verification. The zenith angle changes in the range − 20 to 20 degrees along the square side. Figure 8a presents the matte surface with the distribution of the zenith angle assuming the use of probabilities calculated according to (14) (Lambert law, 2,749,637 registered photons). Figure 8b presents the Gaussian distribution (2,958,408 photons) according to (15) where the dispersion of the zenith angle is set to 10 degrees. Although the Gaussian distribution cannot be verified visually in the same way, the results presented here clearly demonstrate a non-uniform distribution.
The visual verification presented above allowed only qualitative evidence to be obtained. Figure 9 illustrates a comparison between the set distribution function (red line) and the obtained one (points) for the matte (Fig. 9a) and glossy (Fig. 9b) surfaces for the cases presented in Fig. 8. The distributions are obtained along the horizontal lines that pass through the centre. Both obtained distributions are close to set ones with some dispersion influenced by random number generation.

Integration test
The integration test should only demonstrate the correctness and proper action of the entire light propagation model, where the relevant submodels are integrated into one single model. Figure 10 illustrates the results of the integration test. The ratio of glossy and matte reflection is assumed to be 1:1, as well as the length of absorption and scattering, l a = l s . The dispersion of the Gaussian distribution (matte reflection) is set to 5 degrees. Figure 10a visually represents the modelling of light propagation in the integration testing. Figure 10b shows the visual representation of the distribution function when a ratio of matte and glossy reflection is 99:1; a dispersion of the Gaussian distribution is 3 degrees with 1,543,533 registered photons. Figure 10c presents a distribution function along the horizontal line that passes through the centre of the image in Fig. 10b.

Simulation
Heat transfer from the laser beam to opaque material with different reflective properties is analysed and presented in this section.
The laser beam represented by the photons is directed downward on the four spherical particles of the same size, three of which are located in the highest level, and the last is under them. Figure 11a shows a spherical powder particle divided into sectors (there are 1296 sectors in the figure). Four of such powder particles were placed together in such a way that they touched each other, and Fig. 11b presents a snapshot of simulation; photons move downward, struck to the particle surface, then they are reflected or absorbed. A beam of evenly distributed photons strikes the particles from above. The radius of the beam is selected so that it covers all impacted particles in their entirety. When a photon hits a particle, the programme identifies the sector in which the collision took place, and a variable storing the number of photon hits in a given sector of a given particle is incremented. The beam is fired until a given number of hits is reached in any sector (in this case, the number of hits is set to 2,000,000). Photons are reflected by particles, so one photon can hit more than one sector. The number of photon hits in a given sector is divided by the sector area, and as a result, the distribution of the intensity of photon hits on the powder particles is obtained (Fig. 11c).
There are three variants of surface absorption and heat transfer presented in Fig. 11d-f, with one particle hidden for  (Fig. 11e). Most photons reflect multiple times from the surface before absorption in this case. Areas, where probability to accept reflected photon is low, absorb only photons directly falling. The higher the number of different ray tracks leading to the point on the surface, the higher the probability of absorbing photons and the more energy quants are transferred to heat material. The hottest locations are near the particles' contact; some above them, although the contact itself does not receive heat directly from the beam, dark spots can be seen here. The bottom parts of the upper particles are heated from below by reflected light. The absorption and reflection coefficients in the third variant are following: α = 0.05, ρ = 0.95 (Fig. 11f). The view is similar to the previous variant with higher contrast between different areas: upper parts are darker, and locations near contact are brighter. The simulation results can be considered as qualitatively correct.
Described verification methods applied for the developed model confirm adequacy and proper structure and action of the model, which fulfils its functions connected with reflection, refraction, scattering and absorption of laser beam presented by a beam of the photons.

Validation of the model of light propagation
Model validation is based on a comparison of visual and numerical image analysis, theoretical explanation and simulation results. The stainless-steel E-Plus-3D® metal powder was chosen as a representative material for the experiments. The E-Plus 3D 316L stainless-steel powder is a corrosion-resistant iron-based alloy manufactured according to the following standards (ASTM: F3184; DIN EN 10088: X2CrNiMo17-12-2 / 1.4404; UNS: S31603) and applicable for E-Plus 3D's metal 3D printers and other metal printers using the selective power bed fusion technology. As a typical widely known stainless steel with low carbon content that contains in wt% 10-14% Ni, 16-18% Cr, 2-3% Mo and Fe to balance while other elements of its chemical composition were as follows: < 0.03% C, < 1.0% Si, < 2.0 Mn, < 0.045% P, < 0.03% S, < 0.05% O and < 0.1% N. The particle size distribution was within the range of 30-70 μm while other parameters were as follows: angle of repose < 32°; apparent density > 4.3 g/cm 3 and sphericity > 90%. The chosen material has excellent corrosion resistance, heat resistance and creep resistance, widely applied in the food industry, chemical industry, and for mechanical parts, marine equipment, small industrial parts, complex pipes and artefacts.
GX microscope -GT vision, model no. GXM-L3230, with digital camera GXCAM HiChrome HR4 and also CANNON EOS 250D camera, model no. DS126761, has been used for the visual analysis discussed in the following section. Figure 12 presents photographs of L316 steel powder taken using light microscopy at different magnifications. The powder particles are characterised by smooth (indicated by the lines a) and rough (indicated by the lines b) surfaces, as can be seen in Fig. 12a with 1000× magnification. Higher magnification, such as 2500×, allows for observation of more details of particles. The main reflection (Fig. 12b, lines a) is rather large. This may be because of several reasons: the light source is not the point source of light, there is not exact focusing, and/ or the exposition is too high. Peripheral reflection around the main reflection testifies to the high reflection coefficient ρ with the particle surface close to specular or glossy with small dispersion. This is supported by reflections from neighbouring particles indicated by the lines 'b' and by multiple reflections between two particles indicated by the lines 'c' in Fig. 12b.

Visual analysis
At lower magnifications, the powder particles are seen as black spheres with small reflections noticeable as bright points (Fig. 13).
The photographs of the same 316L steel powder taken for different positions of the light source according to the experimental setup are presented in Fig. 14. The first   (Fig. 16). The brightness was calculated as an average value of the brightness of three colours for each pixel. The plots shown in Fig. 16 illustrate changes of the brightness along the two vertical lines from the bottom to the top, the middle line and the line near the right-hand side, i.e. brightness of the powder and brightness of the dark background (matte black paper). The brightness of the powder has low variations without significant trend, especially compared with the background. This allows one to conclude that the reflection is very close to ideal matte (Lambert law).
The reflection of the light photons from the powder was simulated for model verification. A point light source with uniform light emission was located above the powder, similar to the light source shown at Pos. 2 in Fig. 14. The incident photons are reflected without scattering and absorption from the spherical particles of the powder, i.e. simulating single and multiple reflections. The photons directed to the observation area were registered and increased brightness of the appropriate point on the powder surface. Figure 17 presents the simulation results, where an image of the powder layer with a resolution of 400 × 400 pixels is shown in Fig. 17a (with 8,976,052 registered photons) while the distribution of photons vs. the zenith angle obtained in the simulation (blue points) and a theoretical line corresponding to Lambert's emission law (red line) are plotted together in Fig. 17b. There is good agreement, both visual (Fig. 15) and numerical (Fig. 16), between the simulation results obtained and the relevant experimental data.

Theoretical explanation
Matte reflection of light from the layered spherical particles, Figs. 15 and 16, with the specular reflection (Fig. 12) can be explained theoretically. Let us consider the specular reflection from a single spherical particle. A ray of light that falls vertically reflects with the zenith angle, ϑ, which can be described as follows: Then, the distribution function of the reflection will be: Comparing (21) with (5) and (9), it can be concluded that such a reflection is an average of uniform and Lambert reflections in the range of the zenith ϑ angle between 0 and π/2, ( cos 2 = 1+cos 2 ). If the reflected angle is ϑ > π/2, it leads to multiple reflections, which are difficult to analyze analytically, but the same approach can be applied for this reflected light, and it should give a similar result. The result is rather closer to the Oren-Nayar reflection model [46], suggesting that the surface brightness in this case is less dependent on the direction of incident light compared to the Lambert law. Taking this into account, the specular reflection of the wide beam of photons from the single sphere has been simulated for the model verification. The results presented in Fig. 18a illustrate an image of the reflected and 7,056,302 registered photons with a resolution of 400 × 400 pixels while Fig. 18b presents the distribution of the number of photons vs. the zenith angle along the horizontal line crossing the centre of the image. The blue points represent the simulation results, and the red line corresponds to the distribution function calculated according to (21) that shows good agreement.
Fulfilled validation confirms the adequate functioning of the developed model. Further testing of the model requires the implementation of the LBHM into the holistic model of at least one pass of the SLM, the introduction of appropriate model parameters and comparison with the real process.
The first results of model implementation and simulation are presented in the following section, but the final validation is expected later.

Implementation of the LBHM into the SLM model
The model of light propagation discussed in the previous sections has been implemented in the holistic modelling platform of the SLM process characterised by melting and solidification of the deposited powder material and associated energy transfer described briefly in Section 2. The physical events under consideration are powder bed deposition, sin d d laser light propagation and energy absorption, heating of the powder bed by the moving laser beam leading to powder melting, fluid flow in the melted pool, flow through partly or not melted materials and finally to material solidification. The modelling platform can be applied to different multi-material additive manufacturing processes where energy transfer including solid-liquid phase transformation is essential. The code is developed to model two-and three-dimensional cases and deals with powder materials of different particle shapes, sizes and granulometry (distribution). Most of the metal powders used in ALM can be considered atomised (spherical or spheroidal) for modelling purposes, while ceramic or glass powder can be described as granular material, opening new modelling opportunities such as the assumption of different deposition sequences of particles with various shapes and sizes.
The first results presented in Fig. 19 were obtained 4 years ago and have a qualitative character. The simulations were performed on a PC with an Intel Core i7-3930K with the sequential code written in FORTRAN. The modelling space was 150 × 90 cells. A single simulation lasted about 3 h. Figure 19 demonstrates two cases of the simulation of the multi-material SLM process. In both cases, the vertically orientated laser beam moves from the left to the right side of the rectangular domain, transferring energy to the deposited particles. For these examples, bioactive glass and commercially available Ti-6Al-4V powder materials that have Fig. 19 Examples of the singleand multi-material SLM process simulations, a bioactive glass, b bioactive glass and metal. Colours: bioactive glass in a solid state, grey; bioactive glass in a liquid state, blue; metal, black; laser beam, green; intensity of heat transfer, pink-red  Fig. 19, the intensity of light, or the power of the light transferred per unit area, is presented by the intensity of green colour. The heat source, as an absorbed part of the light energy, is shown in pink-red. The solid material is shown as grey for bioactive glass and black for metal, whereas the liquid state is presented as a blue-coloured area. The transparency coefficients of both materials depend on many chemical and physical characteristics of both the laser radiation and the treated materials. Among them are the wavelength of incident radiation, the state of matter, the temperature and also the contributing microstructural defects such as pores and grain boundaries. In the model, the transparency of the particles is varied. They are assumed to be almost opaque in the solid state and much more transparent when the state is changed from solid to liquid. Changes in the transparency coefficient can be adjusted accordingly with available experimental data. The transparency changes near the phase-transformation temperature. The two areas of more intensive heat absorption can be observed in Fig. 19. The first area is located near the liquid surface, where the maximum energy flux density is observed, despite the high Fig. 20 Examples of the multistage and multi-material SLM process simulations, a, b and c the first, second and third layers of Ti-6Al-4V alloy, d the fourth layer of bioactive glass. Colours: unprocessed metal, olive; metal after solidification, brown; liquid, blue, unprocessed bioactive glass, light grey; laser beam, green

Fig. 21
Examples of cross sections at the axis of initial moving laser beam of the third stage of the SLM process simulations (Fig. 20c) showing the temperature distribution: a at the beginning of the stage, b at the end of the laser beam movement to the right, c at the end of the movement to the left, d at the middle of the forward movement, e temperature scale with the real and modelled units. The power intensity of the laser beam is shown by a grade of green transparency of the molten material. The second area is observed near the phase-transition boundary, accompanied by a decrease in material transparency. It can be noticed that the size of the melted pool, its depth and its length depend on the treated material, among other parameters such as the power of the laser beam.
Then, several submodels have been recoded for parallel computations on GPU with CUDA technology using the GPU (initially, it was a GeForce 1060 graphics card with 1280 CUDA cores). The application of such parallel calculations allowed for acceleration at least 100 times. Such simulations of the entire processing cycle do not require complicated interfaces between different components of the model. Some results of the 3D simulation can be found in [38,39], others are presented in Figs. 20 and 21. Figure 20 illustrates the simulation results of the four stages (layers) of the SLM process. The calculation module responsible for simulation of such phenomena as free surface flow, heat conduction and transfer, convection, evaporation, surface tension, wettability with hysteresis and others were developed based on the use of the LBM. Representative volume element contains 256 × 256 × 80 = 5,242,880 nodes (cells) with cell size 1.25 μm, that is, the real sizes are 320 × 320 × 100 μm 3 . The simulation was performed on a PC with CPU AMD Ryzen Threadripper 3960X 24-Core, 4.00 GHz, RAM -32.0 GB, GPU NVIDIA GeForce RTX 3080 and laptop ASUS ROG Strix G533ZS with CPU Intel Core i9-12900H 2.50GHz,RAM 32 GB, GPU NVIDIA GeForce RTX 3080. Parallel calculations were performed using the authors' own code using NVIDIA CUDA, version 11.8 in Microsoft Visual Studio Enterprise 2022, version 17.5.4, and visualisation with the use of OpenGL.
After deposition of the Ti-6Al-4V powder with 1562 atomised particles, three parallel tracks with forward movement of the laser beam were simulated (Fig. 20a). The simulation of one track is approximately 25 min, the cooling after the track is approximately 3 min, all the calculation time is 85 min (1 h 25 min). The unprocessed particles of the Ti-6Al-4V powder are shown in olive colour, the melted material is in blue and the solidified material is in brown. The laser beam is shown in transparent green. Then, the building chamber was filled with the same material (552 new particles). The laser beam moves along the close quadratshaped track in the second stage (Fig. 20b). The powder in the middle of the space remains unprocessed. The calculation time for the second stage is 108 min (1 h 48 min). The third stage (Fig. 20c) is almost the same as the second stage, only a direction of laser movement is opposite. The calculation time of the third stage is 111 min (1 h 51 min). After the third stage, the unprocessed powder was removed from the building chamber and filled with the bioactive glass powder with the 1657 particles of irregular arbitrary shapes. The glass particles are shown in light grey (Fig. 20d). The laser beam does not move in the fourth stage and works in the pulse regime, filling the hole in the middle of the Ti-6Al-4V product with the moulted bioactive glass. The calculation time of the fourth stage is 22 min. In contrast to results presented in Fig. 19, sequential operations with materials with highly different thermal properties give more accurate control over the shape and size of each component of the multi-material product. Figure 21 presents the temperature field in the cross section at the axis of the moving laser beam obtained in simulation of the third stage of the SLM process presented in Fig. 20c. There are four snap shots presented: at the beginning of the stage (Fig. 21a), at the end of the laser beam movement to the right (Fig. 21b), at the end of the movement to the left (Fig. 21c) and at the middle of the forward movement (Fig. 21d).

Conclusive remarks
A laser heating model has been developed as a part of the holistic platform for modelling the multistage multi-material selective laser melting process. The laser heating model is based on the ray-tracing algorithm, which is widely used in visualisation methods and recently for laser heating models in the simulations of selective laser melting. The novelty of the presented model adapted for laser heating in the simulations of selective laser melting is its extension to transparent and translucent materials, including materials with transparency dependent on material temperature. The phenomena of refraction, scattering and volume absorption were implemented into the developed model in addition to the phenomena of reflection and surface absorption commonly considered in such models. The laser beam is considered as a stream of photons of the same energy moving straight without scattering or aberration, whereas reflection, absorption or transmission is randomly chosen on the boundary between different materials with probabilities that depend on the physical characteristics of the material. The direction of reflection or refraction of the photons is defined with the bidirectional scattering distribution function being considered two independent bidirectional reflectance and transmittance functions. The functions defining how light is reflected or transmitted allow for the calculation of the ratio of the light that exits along the outgoing direction to the transmitted light that falls in the incoming direction.
The possibilities and benefits of the proposed approach are demonstrated through numerical test cases representing reflection from a mirror, glossy and matte surfaces and verified by different methods such as units testing, proof of correctness and integration test, while validation of the developed model was based on a comparison of visual and numerical image analysis, theoretical explanation and simulation results. The approach is generic and can be applied to different multi-material SLM processes, where thermal interaction between a laser beam and a powder bed, including solid-liquid phase transformation, is essential. This paper presents examples of simulation of the SLM processes with the developed laser heating model, which demonstrates the possibility of the model in the modelling of complex multistage multi-material processes.
Further development of the platform with the laser heating model is connected with the modification of the SLM LBM-based model in the direction of the SLS/SLM model, with the implementation of a model of the properties of other materials such as bioactive glass, as well as its verification and validation.