Fast Analytic Simulation for Multi-Laser Heating of Sheet Metal in GPU

Interactive multi-beam laser machining simulation is crucial in the context of tool path planning and optimization of laser machining parameters. Current simulation approaches for heat transfer analysis (1) rely on numerical Finite Element methods (or any of its variants), non-suitable for interactive applications; and (2) require the multiple laser beams to be completely synchronized in trajectories, parameters and time frames. To overcome this limitation, this manuscript presents an algorithm for interactive simulation of the transient temperature field on the sheet metal. Contrary to standard numerical methods, our algorithm is based on an analytic solution in the frequency domain, allowing arbitrary time/space discretizations without loss of precision and non-monotonic retrieval of the temperature history. In addition, the method allows complete asynchronous laser beams with independent trajectories, parameters and time frames. Our implementation in a GPU device allows simulations at interactive rates even for a large amount of simultaneous laser beams. The presented method is already integrated into an interactive simulation environment for sheet cutting. Ongoing work addresses thermal stress coupling and laser ablation.


Introduction
Multi-beam laser heating of sheet metal is a relevant metalworking technique which has arisen interest of researchers in the last years. In contrast to single-beam heating, multi-beam heating provides two main advantages to the former: (1) the ability to process different locations of the sheet simultaneously [1], and (2) control of thermal stress levels by specific multi-beam configurations [2]. Industrial applications of multi-beam heating of sheet metal include laser forming and bending, laser cutting and additive manufacturing.
Thermal simulation is crucial for temperature and stress analysis of manufactured pieces. An adequate selection of laser parameters and a correct path planning allows for improving the efficiency of the process and minimizes material damage and waste. Current simulation approaches rely on numerical schemes which require fine geometry and time discretizations. Such discretizations imply high computational costs, which limit the application of these approaches in real manufacturing scenarios with large time/space domains and complex laser trajectories [3].
This manuscript presents a simulation approach for the multi-beam laser heating problem based on an analytic solution to the heat equation on rectangular domains. This analytic solution does not require a mesh nor a fine time discretization to solve the problem. As a consequence, our algorithm is able to run complex simulations with large time/space domains and complex multi-laser trajectories at interactive time rates. Furthermore, each laser beam trajectory is allowed to be independent from the others, with different time-dependent parameters, trajectories and time frames (i.e., each laser beam can be turned on/off independently at any point of the simulation).
The remainder of this manuscript is organized as follows: Section 2 discusses the relevant literature. Section 3 presents the mathematical model and describes the implementation of the proposed algorithm. Section 4 discusses the obtained results for different test cases. Finally, Section 5 presents the conclusions and introduces the future work.

Multi-Beam Single Trajectory vs. Multiple-Trajectory Simultaneous Laser Heating
There are currently two main applications for multi-laser heating in laser machining: (1) single-trajectory multi-laser heating and (2) multi-trajectory multi-laser heating.
In single-trajectory multi-beam laser heating, a leading laser is followed by a pattern of secondary finishing ones, all in the same trajectory (or with minimum spatial offset). Experimental evidence has shown that such a configuration reduces the thermal stresses produced by the laser beams in laser cutting (compared to single-beam cutting) [4]. Furthermore, specific configurations of the laser beams have been shown to reduce the required pressure of assisting gas [5]. Each laser beam may be produced by an independent source [6] or by diffraction of a single beam source [7].
On the other hand, in multi-trajectory heating, each laser beam follows an independent trajectory [1]. Multi-trajectory heating is relevant as it improves the machining times by processing different zones of the sheet at the same time. This manuscript focuses on the simulation of multi-trajectory laser heating.

Thermal Simulation for Sheet Metal Laser Heating
The problem of laser beam heating simulation has been widely researched for single-beam applications. Numerical methods are the most common simulation approach. Methods such as the Finite Differences Method (FDM) [8][9][10][11] and the Boundary Element Method (BEM) [12,13] have been used in the literature to study the thermal behaviour of sheet laser heating. However, these methods impose several numerical limitations such as dense rectangular grids in the case of FDM and non-sparse linear systems for BEM, requiring a high amount of computational resources even for simple problems.
The Finite Element Method (FEM) is a standard numerical approach for the simulation of physical problems. Nonlinear FEM has been applied to study the thermal/stress behaviour of the single-beam laser heating of rectangular sheets [14][15][16][17][18]. In contrast to FDM and BEM, FEM works by discretizing the domain using different types of meshes which allow fine discretizations (high level of detail) near interest zones (laser spot and trajectory) and coarse resolutions in other zones resulting in less expensive systems of equations. To address laser ablation and material removal, methods such as element birth and death [19], volume fractions [20], temperature thresholds [21][22][23], and the enthalpy method [15][16][17][18]24] are coupled to the different numerical schemes.
The Finite Volume Method (FVM) has been recently introduced in the literature for the study of laser ablation and sheet heating [25,26]. Contrary to previous numerical methods, the FVM allows to accurately model and simulate interactions between the laser beam and the sheet in its physical states: solid, liquid and gas. However, these interactions are highly nonlinear and, as a consequence, computationally expensive.
Numerical methods provide tools to simulate nonlinear interactions, phase changes, laser ablation and material removal. However, they are also computationally expensive, limiting the applications of the algorithms in real manufacturing scenarios with complex laser trajectories and large space/time domains. Analytical methods do not have such limitations, allowing fast simulation of complex problems at the cost of some simplifications of the physical model. Uni-dimensional analytic models have been implemented for the simulation of laser drilling processes for static laser beams [27,28].
A 3D model for laser heating of sheet metal for straight line trajectories is presented in [29], while 2D models for arbitrary laser trajectories have been presented in [30,31].
Despite the large amount of literature concerning single-beam laser heating, simulation for multi-beam laser processing has been rarely studied. In [6,32], the FDM is used to analyze the thermal and structural impact of two independent laser beams melting a sheet metal. On the other hand, in [2,33], a thermal/stress analysis of multi-beam laser heating is performed using FEM. These simulations show that multi-beam heating reduce the thermal stresses afflicted to the material. A semi-analytic model for the steady multi-beam laser heating problem is presented in [34]. This pseudo-analytic model is also implemented inside an optimization algorithm which estimates the best laser parameters for a given manufacturing process. The use of analytic solutions is crucial for optimization due to the optimization process being expensive per se, requiring multiple evaluations of the temperature fields with different laser parameters and/or trajectories.

Conclusions of the Literature Review
Contrary to single-beam, multi-beam heating is scarce in the literature. Numerical (as opposed to analytical) methods are computationally expensive. As a consequence, these methods are unable to simulate real world sheet sizes and laser trajectories.
This manuscript offers the implementation of an analytic Fourier-based method to simulate multi-trajectory laser heating. The characteristics of the implemented method are: (a) constant material properties, (b) natural convection, (c) simplification of the sheet into a 2D domain, (d) transient (time-history) temperature, (e) multiple laser head configurations, (f) independent laser trajectories (with independent parameters), and (g) arbitrary power on/off time history on each trajectory.
Our algorithm is heavily based on the previous work presented in [31]. However, our method presents several improvements to this previous work as follows: (a) our algorithm reformulates the analytic solution to allow simultaneously multiple laser beams instead of only one. (b) Our method allows completely asynchronous laser beams (independent parameters, trajectories and time frames) which is not a problem for single-beam approaches and has not been addressed in other multi-beam related work. (c) Instead of computing and storing the full temperature history, our implementation exploits the frequency-based solution approach by interactively providing the temperature only when requested by the user (even in non-monotonic order). (d) Computing the solution in the frequency domain allows arbitrary sheet discretizations for visualization and analysis, including triangular meshes, rectangular grids, polylines or sets of sampled points on the sheet. (e) Our GPU implementation of the algorithm using on-chip (local) memory presents an increased 10x speed-up on computation time over previous work.

Heat Transfer Equation for Multi-Beam Laser Heating
The temperature distribution u = u( x, t) of a 2D rectangular sheet satisfies the following partial differential equation: where ρ, c p and κ are the material density, specific heat and thermal conductivity, respectively. ∇ · ∇ = ∂ 2 /∂ 2 x + ∂ 2 /∂ 2 y is the 2D Laplace operator. ∆z is the thickness of the sheet. F = F( x, t) is the set of surface heat sources affecting the sheet, and q = q(u) is the temperature-dependent convection on the sheet surface with its respective convection coefficient h ≥ 0 and ambient temperature u ∞ ∈ R. For the previous partial differential equation, the following boundary and initial conditions are imposed on the sheet: For the multi-beam approach, the set of heat sources F is defined as the sum of the heat inputs of each laser beam f k : where P k ≥ 0 is the power of the laser beam, r k > 0 is the radius of the laser spot and x k 0 (t) ∈ R 2 is the center of the laser spot at time t. 0 ≤ R < 1 is the reflectivity of the material. x ∞ = max(x, y) is the infinity norm in R 2 . The above laser model transforms the circular laser spot with radius r k to its equivalent squared spot with area πr 2 k . We apply such transformation in order to develop the analytic solution in the next section. Figure 1 presents a scheme of the multi-beam laser heating problem.

Laser beams
Convective heat loss

Analytic Solution
Following the same procedure as in [31], the temperature u on the rectangular sheet can be expressed as a linear combination of Fourier functions: with basis: Applying separation of variables, the Fourier coefficients Θ ij from Equation (4) can be expressed as a sum of the pseudo-coefficients θ k ij of each independent heat source: where each laser beam f k defines its respective pseudo-coefficient θ k ij as follows: with corresponding Laplacian eigenvalues: Equation (7) is written recursively in terms of a previous known solution θ k ij (t 0 ). Therefore, each laser trajectory is discretized as a piecewise linear trajectory . In contrast to previous work [31], Equations (3), (6) and (7) have been reformulated to account for the multiple asynchronous laser beams. Finally, the closed form for the integral term in Equation (7) has been already presented in [31]. 1.

Algorithm
Discretize laser trajectories: As discussed in Section 3.2, the laser beam trajectories , as described in [31]. The only requirement for this discretization is the fact that all laser beam trajectories must share the same time discretization, i.e., t 0 , t 1 , . . . , T f are the same for all trajectories x k 0 .

2.
Compute the Laplacian eigenvalues: The Laplacian eigenvalues of the sheet are computed as per Equation (8). Since the eigenvalues ω ij are time-independent, this step is performed before the simulation loop starts.

3.
Initialize time and sheet temperature: The simulation time is initialized to t 0 ← 0. In order to satisfy the initial temperature condition u(0) = u ∞ , the pseudo coefficients are initialized as θ k ij (0) = 0 (see Equation (4)).

4.
Update current time t: The current simulation time t = t l+1 is updated according to the previous time t 0 = t l , in concordance with the discretization of trajectories from step 1.

6.
Question: Is laser beam k turned on?: This step allows for simulating asynchronous laser beams by asking at the current time t if the laser is turned on/off. Therefore, each laser beam might have its own internal time Set power P k /Set null power P k ← 0: In the previous step, the program checks the state (on/off) of the current laser beam k. The laser is turned on by the simulation by setting its corresponding power input P k . In the case of the laser being turned off, the simulation simply sets its power to 0. 8.
Compute the pseudo-coefficients for each laser beam: The pseudo-coefficients in Equation (7) are solved analytically for each laser beam as described in [31]. The number of coefficients computed is truncated to num_coe f f s.

9.
Question: k < num_lasers?: Check if the pseudo-coefficients have been computed for all the laser beams. 10.
Compute the Fourier coefficients as per Equation (6): This step computes the real Fourier coefficents Θ ij (t) from the pseudo-coefficients θ k ij (t) as per Equation (6).

11.
If required, compute temperature: The temperature field is computed on a set of discrete points sampled on the sheet [(x 0 , y 0 ), (x 1 , y 1 ) . . . ] as per Equation (4). The number of coefficients used to compute the solution is truncated to num_coe f f s. This step is optional since the result may be stored in the frequency domain (Θ ij (t)). Therefore, the temperature is made available only when requested by the user, allowing for skipping iterations of no interest and allowing non-monotonic access to the transient temperature history, improving the performance in the process. 12.
t  Figure 2. Diagram of the analytic multi-laser simulation algorithm.
There are several concepts in the previous algorithm that help to improve its efficiency, crucial for the interactive nature of the simulation: 1.
In step 1 of the previous algorithm, the curved laser beam trajectories x k 0 (t) are discretized as piecewise linear ones [ x k 0 (0), x k 0 (t 1 ), . . . , x k 0 (T f )], which inherently produces the time discretization [0, t 1 , . . . , T f ]. This time discretization does not affect the numerical accuracy of the temperature solution. Therefore, as opposed to FEA (Finite Element Analysis), the time step size ∆t = t l+1 − t l of our algorithm can be arbitrarily large. 2.
Step 6 allows turning on/off any laser beam at any point of the simulation, allowing complete asynchronicity between the multiple laser beams. In addition, the algorithm allows to change any laser parameters at will, resulting in time-dependent parameters P k (t) (laser power) and r k (t) (spot radius). For simplicity, this manuscript uses constant parameters. 3.
The complete information of the solution is stored in the frequency domain (step 10) and temperature data is computed only when requested. Therefore, the user requests the temperature only at specific times and in specific zones (i.e., at the middle or the end of the simulation). This frequency-based approach has the advantage of providing non-monotonic access to the temperature history, allowing arbitrary simulation time-location requests. Furthermore, since the space discretization does not affect the solution, any sheet sampling can be used to render the temperature (rectangular grid, triangular mesh, a curve or a single point in the sheet). 4.
In step 10, each pseudo-coefficient θ k ij is independent from the rest of the pseudo-coefficients (Equation (6)). Similarly, in step 11, the temperature value u at a given point x is independent from the temperature on the rest of the sheet (Equation (4)). Therefore, the computation in both of these steps is parallelized.

Results
This section presents the simulation results of our algorithm. For all the simulations, a Fourier discretization of 512 × 512 coefficients and a spatial (grid) discretization of 512 × 512 points are used. These resolutions have been chosen as they have shown in our experiments enough accuracy within our desired execution time ranges (see Sections 4.1 and 4.3, respectively). As the ratio max(a,b) min k (r k ) (sheet size vs. spot size) increases, it could be necessary to increase the number of Fourier coefficients. However, a sensitivity analysis of our algorithm's discretization is out of the scope of this manuscript. Table 1 presents the geometric and physical parameters of the sheet used in the simulations. All of the experiments presented in this manuscript do not consider the effect of heat reflected by the sheet (R = 0). The first study case presents two different laser beams heating the sheet simultaneously. The first laser beam follows a star-shaped trajectory on the sheet while the second laser beam follows a rectangular trajectory. Table 2 presents the parameters of each laser beam. As discussed in Section 3, our algorithm not only enables different parameters for multiple lasers (path, power, speed, spot radius), but it also allows independent time frames for each trajectory. Parameter values for Tables 1 and 2 are taken from state-of-the-art literature [14,31,35]. However, laser beam power has been down-scaled to account for sheet heating instead of cutting or bending. Figure 3 plots the simulation results for two laser beams heating the sheet surface. Figure 3a shows the star and square laser trajectories planned on the sheet. The laser beam parameters for each trajectory are described in Table 2. At the beginning of the simulation, a unique laser beam (star) heats the surface (Figure 3b). As discussed in Section 3.3, our algorithm allows independent time frames for the multiple laser beams. Therefore, the second laser is introduced in the simulation at t = 0.065 s (Figure 3c). Figure 3d plots the temperature when the two lasers are near each other (the trajectories do not intersect). Finally, both laser beam trajectories end at different steps: t k f = 0.13 s for the square trajectory ( Figure 3e) and t k f = 0.16 s for the star trajectory (Figure 3f).

Comparison with FEA
To validate the analytic approach, FEA simulation is performed. The FEA linear system for Equation (1) becomes [31]: where U and F k are the vectors of temperature and heat sources sampled on the sheet, and L and M are the stiffness and mass matrices, respectively [31]. The software ANSYS R Academic Research Mechanical, Release 17.2, is used to perform the FEA simulations. ANSYS R element SHELL131 is employed. Element thickness and material properties are set as per Table 1. The elements are configured to have a constant temperature along the thickness and to evaluate convection at the sheet surface, as per Equation (1). To represent the area heated by the laser beams at every time step, ANSYS R surface loads (Heat Fluxes) are applied on the FEA elements that lie inside the heated zone. Figure 4 plots the FEA results for the study case presented in Table 2. The mesh of the domain is generated so that it is more dense in the neighborhoods of the laser trajectories. Figure 4a shows the initial triangular mesh computed using the FEA software which is then refined several times (Figure 4b) to improve the numerical accuracy of the solution. After five re-meshing iterations, the final mesh contains 126 k triangles and 63 k nodes. Figure 4c plots the FEA temperature at the moment the two laser beams get closer to each other (t = 0.094 s). The temperature peak is in the middle of the two trajectories (Figure 4d), due to both paths not intersecting each other.    (Figure 4c). The maximum relative error is 3.9%, located around the laser spots (Figure 5b). This small deviation between our method and FEA arises because: (a) the number of Fourier coefficients used in implementation of Equation (4) affect our algorithm precision, (b) the precision of the FEA algorithm is dependent on mesh and time resolution, element type, etc., and (c) the FEA approach approximates the laser spots using FEA elements at each time step, not preserving the exact geometry of the laser spots. The small square shape of the error in Figure 5b is due to the squared laser model presented in Equation (3).

Multiple Laser Beams
The algorithm introduced in this manuscript allows multiple laser beams heating the surface at the same time. This section presents additional experiments with more than just two laser beams. Figure 6 presents a study case with four simultaneous laser beams drawing different shapes on the sheet: a square, a star, a spiral and a circle trajectory (Figure 6a). All of the laser beams have the same parameters (power and radius) and start at the same time (Figure 6b). However, they do not finish at the same time (Figure 6c). As discussed in Section 3.3, our algorithm allows for defining different independent time frames for each laser beam by turning on/off specific laser beams. Such approach even allows to turn off all the laser beams and continue the simulation. Figure 6d illustrates this approach by continuing the simulation after all the laser beams have finished their trajectories, where only thermal conduction and thermal convection are taken into account.
In order to visually compare our algorithm with other simulation approaches in the literature, the study case presented in [33] is replicated in this manuscript (Figure 7). In this study case, seven simultaneous lasers are distributed uniformly on the y-axis. Each laser beam follows a straight line trajectory and together they draw an arc in the heating front (Figure 7a,b). Our algorithm is able to produce similar results to [33] despite the simplification of the analytic model (Figure 7c,d).  Figure 8 plots the temperature distribution along the arc for different number of laser beams. As the number of laser beams increase (and the arc length remains the same), the oscillations of the temperature are mitigated and the temperature increases. Such a result is consistent with the analysis presented in [33]. As discussed in Section 3.3, our algorithm allows for computing easily the temperature along the arc without even requiring to calculate the temperature on the rest of the sheet, which improves the computation time for this particular case.  Figure 7. Our results are similar to those in [33]. ns represents the number of lasers.

Performance Assessment
This section evaluates the performance of our algorithm under hardware-accelerated (GPU) and non-accelerated (CPU) platforms. The presented algorithm has been implemented using the OpenCL framework to create an optimized solution that targets both CPU and GPU. On systems where only a CPU is available, our implementation makes use of multi-core parallelization and vectorization to speed-up computation. On systems where a GPU is available, the memory hierarchy can be explicitly controlled using the OpenCL API. The workload is divided into small groups, in order to exploit reuse of computed data using local memory. In this manner, a high speed-up is achieved due to both efficient use of memory and massive parallelization.
The performance results have been measured with the following test platform: A desktop PC using Windows 10 OS with an Intel R Core TM i5-6500 (CPU), 8 GB RAM and NVIDIA R GeForce R GTX 960 graphics card. Our algorithm is able to simulate any number of laser beams. Figure 9 plots the execution times for the computation of the Fourier coefficients as a function of the number of laser beams. The figure compares the execution times of the CPU against the GPU to compute 512 × 512 Fourier coefficients. The computational cost increases with a large slope in the CPU approach while being nearly constant in the GPU approach. The more laser beams are added, the more it benefits from GPU parallelization. In addition, the computation of the Fourier coefficients in the GPU is preferable. In this way, there is no need to transfer the coefficients back and forth from host-to-device on each iteration since they always stay in GPU memory. For this analysis, a single time step has been considered instead of the whole simulation.  Figure 10 shows the computation times of the temperature retrieval in a mesh grid as a function of the grid size. While the curves for Fourier resolutions (number of coeffs) below 512 × 512 display a computational cost relatively low (<0.5 s) for any spatial resolution, Fourier resolutions above that point (1024 × 1024) become expensive (>0.5 s) for our simulation purposes. Therefore, on our test platform, we have observed that a good balance between accuracy and computation time can be achieved by setting a resolution of 512 × 512 for both the frequency (Fourier) and spatial (grid) domains. Compared to the previous work [31], a 10x speed-up is achieved at the temperature evaluation by using local GPU memory. Additionally, since the computation of the temperature is not compulsory at each iteration of the algorithm (as discussed in Section 3.3), the simulation may ignore irrelevant time steps (defined by the user). Moreover, the independence of the algorithm accuracy from the spatial discretization allows the user to specify specific domains (e.g., a sheet section, a curve or a single point) without requiring the whole sheet temperature. These aspects improve the computation efficiency of the simulation for specific requirements of the user.

Integration within an Interactive Laser Cutting Simulation Environment
Nowadays, multi-laser machines with multi-trajectory capabilities do not represent a significant share of the market. The current state-of-the-art laser cutting machines can be divided in three main groups: (1) multi-laser machines with independent and disjoint working areas for each laser head, (2) multi-laser machines with parallel heads working simultaneously and (3) multi-laser machines with fully individual and autonomous laser heads (whose trajectories may intersect or get close enough). The first scenario can be reduced to a collection of mono-laser machining scenarios, since the individual working areas for each laser head avoid interference with the neighbouring laser heads and their heat effects. Therefore, the approach introduced in [31] is still valid to simulate the temperature for each laser head of these multi-laser machines.
In the second type of multi-laser machines, all laser heads receive the same trajectories and machining commands, but their separation or offset can be setup and reconfigured during the machining. Additionally, each laser head can be switched off and on individually. The offset between the heads does not guarantee disjoint working zones. Therefore, the simulation approach presented in this work is used to address the potential interference between the multiple heat sources.
The last type of multi-laser machines use mirrors and lenses to move the spot of the available laser heads. Each laser spot can be commanded independently of the others, even allowing two or more spots to converge at the same physical point on the sheet metal. The methodology presented in this work is used to simulate such scenarios.
In the context of the second scenario, this section presents an interactive simulator of a laser cutting machine with three laser heads. All laser heads are arranged side by side in the bridge of the machine (x-axis of the machine). Each individual laser head can be switched on and off individually and their relative positions among them (y-axis of the machine) can be changed by means of specific machining instructions.
The virtual simulator itself uses a contour-based representation to model the geometry of the processed sheet [36,37], presenting a virtual 3D interactive scenario with the multi-laser CNC (Computer Numerical Control) machine that receives the machining instructions. In this virtual scene, the moving and cutting instructions move the corresponding axis of the machines (bridge, laser heads offset along the bridge and laser head height over the sheet metal). The heat sources are then calculated and sent to the heat simulation procedure, updating the temperature of the sheet, which is rendered as a texture over the virtual sheet metal.
Since all the laser heads receive the same machining instructions (although their position along the bridge differ), all movements start and end at the same time. As discussed in Section 3.3, the laser beam trajectories must share the same time discretization, which in this case, is guaranteed by design. If a unrestricted multi-laser machine is used with independent laser heads, i.e., receiving different machining instructions, a global time clock must be used in order to trigger the update of the temperature computation with the correct positions of the moving laser heads.
The introduction of the multi-laser machines in the industrial world has an impact from the design point of view. Even with just one laser, the designer of the NC programs must take into account the expected produced heat and how it spreads over the sheet metal in order to optimize the nesting of the produced parts. With a multi-laser machine, this procedure is even more critical, as the resulting heat of the multiple sources can accumulate in some areas. NC designer is expected to use the presented virtual simulator to visualize how the cutting process produces the parts and to analyze the computed temperature through the sheet. If any situation becomes problematic, the designer would make changes to the NC program in order to address the situation.
The presented multi-laser simulator runs at interactive rates. Therefore, the designer can improve the optimization workflow, as many simulations can be run in a short period of time. In contrast, high quality simulations with FEA software, although being highly precise, are computationally expensive. Thus, the number of test configurations that the designer can test during the design phase is limited. Nevertheless, at the end of the optimization phase, the interactive heat simulations can be complemented with high quality FEA simulations.
From the industrial point of view, this improved design workflow, assisted with the interactive multi-laser simulation, (1) provides better NC programs in terms of the quality of the produced parts, (2) produces economical benefits due to shorter machining times or less wasted resources, and (3) improves the safety of the operators in the factory floor, reducing unnecessary risks due to unexpected behaviour of the NC programs. Figure 11a shows a machine with three laser heads working simultaneously, i.e., they receive the same machining instructions with fixed spatial offsets). Each laser head machines a star figure. The stars machined by the first and second laser heads overlap, while the third laser head produce an isolated star figure. Figure 11b shows a closer view of the cutting area while Figure 11c shows the same viewpoint, but the cut stars removed from the visualization. The simulation results show that intersecting trajectories present temperature peaks (black zones) where the trajectories overlap (see Figure 11b).

Conclusions
This manuscript presents a novel methodology for the simulation of heat transfer on rectangular sheet metal under multi-laser beam heating. The algorithm is based on an analytic frequency-based solution to the heat transfer equation which considers some simplifications (2D domain, constant material parameters) in favor of simulation speed. Such simplifications allow the algorithm to solve the transient temperature map on large space/time domains and complex laser trajectories. Furthermore, the algorithm allows many simultaneous laser beams with independent parameters (laser power, speed and spot radius) and time frames (i.e., each laser beam can be turned on/off during the simulation). The numerical accuracy of the algorithm is independent from the space/time discretization, allowing non-monotonic access to the temperature time history and arbitrary discretization of the sheet domain. Our simulation tests show that the algorithm is able to render correctly the temperature maps for several laser beams with different trajectories using mesh grids.
A numerical comparison with FEA shows that our algorithm solution deviates from the FEA one a maximum of 3.9% and only around the laser spot. An assessment of the algorithm performance shows that, in an implementation using GPUs, the number of laser beams barely affects the simulation time. In addition, compared to previous work [31], the temperature retrieval time has been reduced by a factor of ten. Finally, the algorithm is implemented into an interactive laser cutting simulation environment for the assessment in real time of laser cutting processes in real manufacturing scenarios.
The presented algorithm simplifies by design the mathematical model of the problem in favor of interactive simulations. As future work, we work on: (1) simulating the nonlinear behaviour of the material properties which arises due to the high temperature gradients, (2) simulating physically the laser ablation and material removal in sheet cutting processes, (3) coupling the model with an analytic stress model in order to evaluate the potential structural damage due to thermal stress, and (4) studying the sensitivity of the Fourier discretization with respect to parameters such as sheet size and radii of the laser spots.

Acknowledgments:
The authors thank Lantek Business Solutions and Lantek Investigación y Desarrollo for providing the laser cutting machining NC programs used in this work.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: Temperature distribution along the sheet at any given time (K).
κ Sheet thermal conductivity (W/(m K)). R Sheet reflectivity (0 ≤ R < 1). q = q(u) Temperature-dependent heat convection on the sheet surface (W/m 2 ). u ∞ Ambient temperature K. h Natural convection coefficient (W/(m 2 K)). f k Heat input from laser beam k (W/m 2 ). k = 1 . . . num_lasers. P k Power of laser beam k (W). r k Radius of laser spot k (m). x k 0 = x k 0 (t) Laser spot center for laser beam k at time t. [t k 0 , t k f ] Simulation time frame in which the laser beam k remains turned on. 0 ≤ t k 0 < t k f ≤ T f . v k Scan speed of laser beam k (m/s).
Sum of all laser beam heat sources (W/m 2 ).
X i = X i (x) Fourier basis function associated to the x coordinate. i = 1 . . . ∞.
Y j = Y j (y) Fourier basis function associated to the y coordinate. j = 1 . . . ∞. Θ ij = Θ ij (t) Fourier coefficient associated to basis functions X i and Y j . θ k ij = θ k ij (t) Pseudo-Fourier-coefficient associated to the k-th laser beam source.
ω ij ij-th eigenvalue of the heat operator (Laplacian) on the rectangular sheet.