Efﬁcient Computation of Heat Distribution of Processed Materials under Laser Irradiation

: In this paper, a solution is provided to solve the heat conduction equation in the three-dimensional cylinder region, where the laser intensity of the material irradiation surface is expressed as a Gaussian distribution. Based on the symmetry of heat distribution, ﬁrstly, the form of the heat equation in the common rectangular coordinate system is changed to another form in the two-dimensional cylindrical coordinate system. Secondly, the ADI with the backward Euler method and with Crank–Nicolson method are established to discretize the model in the cylindrical coordinate system, after which the simulation results are obtained, where the ﬁrst kind of boundary value condition is used to verify the accuracy of these two algorithms. Then, the above two methods are used to solve the model with the third kind of boundary value condition. Finally, the comparison is performed with the results obtained by the MATLAB’s PDETOOL, which shows that the solution is more feasible and efﬁcient.


Introduction
In the fields of advanced equipment manufacturing for example aerospace and new energy, hard and brittle materials such as beryllium, fused silica, and diamond are widely used to manufacture products and devices. Generally, when machining (drilling or cutting) brittle and hard materials, it is pretty easy to cause damage on the processed material. Therefore, at present, laser is commonly used to irradiate the surface of the processed material, after which the processed materials obtain heat through the interaction between light and itself, and the processing quality is significantly improved by heating and modification. In order to find the optimal laser intensity and distribution, it is necessary to calculate the propagation state of light and heat in three-dimensional object. Common numerical methods are the finite element method (FEM) and the finite difference method (FDM) [1,2], but they all require large amounts of computation.
For the purpose, some numerical methods had been proposed [3][4][5][6][7][8][9][10][11][12][13][14], such as the alternating direction implicit (ADI) finite-difference time-domain method (FDTD) and its convolution perfect matching layer (CPML). The implementation of the ADI-FDTD and its CPML is divided into three steps. Firstly, transform the finite difference time domain method in the traditional three-dimensional cylindrical coordinate system into a matrix expression. Secondly, matrix expression of ADI-FDTD method in three-dimensional cylindrical coordinate system is proposed by matrix transformation. Finally, the influence parameters of the marching layer are added to the method in the form of auxiliary variables.
The existing methods directly give a three-dimensional cylindrical coordinate model in the electromagnetic environment, combined with the unconditionally stable ADI algorithm. However, we propose a three-dimensional heat conduction model of laser processing in the rectangular coordinate system, which is then discretely transformed, simplified, and solved in three steps. Firstly, a three-dimensional heat conduction model in a rectangular coordinate system is developed, and it is transformed into a three-dimensional heat conduction model in the cylindrical coordinate system. Secondly, since the heat conduction discussed in this article is irrelevant with the angle in the cylindrical coordinate system, the threedimensional problem in cylindrical coordinates is simplified to a two-dimensional problem. Finally, we use the backward Euler method and Crank-Nicolson method, combined with the unconditionally stable ADI method for discretization.
In order to solve this problem, we first develop the heat conduction equation and boundary conditions. According to the characteristics of laser beam, the heat conduction equation solving problem in three-dimensional rectangular coordinate system is transformed into the heat conduction equation solving problem in two-dimensional cylindrical coordinate system by using cylindrical coordinate transformation. Finally, an alternate implicit scheme algorithm is constructed. In conclusion, the heat conduction distribution can be solved quickly and stably, which provides an effective calculation method for the optimization of related parameter.

Mathematical Model
The problem model discussed in this article is first proposed in a rectangular coordinate system, and then converted to a corresponding cylindrical coordinate system for discretization and solution. The significance of this chapter is to give the origin of the cylindrical coordinate equation model in the article and connect it to the original equation model in the rectangular coordinate system.
Suppose Ω is a cylinder with the origin of the coordinates as the center, then the model of the problem for three-dimensional heat conduction is as follows: T(x, y, z, 0) = ϕ(x, y, z), 0 ≤ x ≤ X, where α ≥ 0. Boundary value conditions of the first kind are: T S 1 = k(y, z, t), S 1 = (x, y, z)|x = 0, y 2 + z 2 ≤ R 2 , T S 2 = l(x, t), S 2 = (x, y, z)|0 < x < X, y 2 + z 2 = R 2 , T S 3 = g(y, z, t), S 3 = (x, y, z)|x = X, y 2 + z 2 ≤ R 2 . (3) Boundary value conditions of the third kind are: κ ∇T·n Γ = βq y 2 + z 2 , Γ = (y, z)|x = 0, y 2 + z 2 ≤ d , where Γ is the laser irradiation area, d is the laser radius, β is the material's absorption rate of laser energy, H is the heat dissipation coefficient, and T a is the initial environmental temperature, as shown in Figure 1. Consider the problem model under the column coordinates and make the following coordinate transformation: Then, Equation (1) is transformed into the following form: whereT(x, y, z, t) = T(x, r, t).
When the temperature distribution T is independent of the θ, Equation (6) is transformed into the simpler form: Transforming Equation (1) into Equation (7) is equivalent to transforming threedimensional problems into two-dimensional problems. Compared with the direct solution of three-dimensional problems, we have simplified the problem, saved time, and improved the efficiency.
Consider the transformed two-dimensional problem model as follows: Boundary value conditions of the first kind are as follows: Boundary value conditions of the third kind are as follows: The alternate direction implicit format (abbreviated as ADI format) is an unconditionally stable format which can be solved by the catch-up method. This method is also called the P-R method, which was proposed by Peaceman and Rachford in 1955, and will be used to calculate the matrix form under the column coordinate.

Establishment of Differential Approximation
In the case of X = R, the space step is set to be h = (X − 0)/M, and the time step is set to be τ = (T e − 0)/N. Where M and N are positive integers, we have x i = 0 + i·h, 0 ≤ i ≤ M, r j = 0 + j·h, 0 < j ≤ M, and t k = 0 + k·τ, 1 ≤ k ≤ N. So we get a grid subdivision of the interval, and x i , r j , t k is a node.

Difference Scheme for the Backward Euler Method
At node x i , r j , t k+1 , Equation (8) is considered as: Using the backward Euler method to discretize Equation (12), we obtain [2]: where the total truncation error is as follow: The difference equation can be obtained by replacing the exact solution T x i , r j , t k+1 with an approximate numerical solution T k+1 i,j and discarding the truncation error R (1) i,j,k :

Difference Scheme for Crank-Nicolson Method
Let f = 0, and at node x i , r j , t k+1/2 , homogeneous form of Equation (8) is considered as the following form: Using Crank-Nicolson method to discretize Equation (16), we obtain [2]: where the total truncation error is: The difference equation is obtained by replacing the exact solution T x i , r j , t k+1 with a numerical solution T k+1 i,j and discarding the truncation error R (2) i,j,k : The difference scheme can be obtained by synthesizing the initial value condition and the boundary value condition as follows: Denote: The difference equation is written as: There exists a normal number C 1 independent of h and τ that causes the total truncation error: In Equation (23), adding the minor term α 2 τ 2 δ 2 x δ 2 r + 1 r j δ r T k+1 i,j − T k i,j leads to: To facilitate the calculation, the transition layer variable T k+1/2 i,j is introduced, and the above equation becomes: Denoter = α·τ h 2 andr 2 = α·τ 2h , then the vector form of the above two equations is: Denote Denote b 1 (j) =r +ˆr 2 r j = 1 + 1 2j ·r, b 2 (j) =r −ˆr 2 r j = 1 − 1 2j ·r, then we have: The transition layer variable T k+1/2 0 , T k+1/2 M should satisfy: There exists a normal number C 2 , independent of h and τ, so that the minor term is: Therefore, Equations (24) and (38), show that a normal number C is independent of h and τ, so the total error R of difference schemes for ADI with the backward Euler method of the first kind of boundary value is:

Difference Scheme for ADI with Crank-Nicolson Method
The difference scheme is obtained by synthesizing the initial value condition and the boundary value condition: Denote: The difference equation is written as: There exists a normal number C 1 independent of h and τ that causes the total truncation error: In Equation (43), adding the minor term α 2 · τ 2 4 ·δ 2 To facilitate the calculation, the transition layer variable T k+1/2 i,j is introduced, and the above equation becomes: Denoter = α·τ h 2 , then the vector forms of the above two equations are as follows: then we derive following forms: ·r, then we get: The transition layer variable T k+1/2 0 , T k+1/2 M should satisfy: There exists a normal number C 2 , independent of h and τ, so the minor term is: Therefore, Equations (44) and (56), demonstrate that the normal number C is independent of h and τ, so the total error R (1) of difference schemes for ADI with Crank-Nicolson method of the first kind of boundary value is: Based on the initial value condition, the boundary value condition and Equation (25), the difference scheme is obtained: It is the same as the difference scheme under the boundary values of the first kind, and finally discretized into the corresponding ADI scheme. Therefore, the rest are consistent, except that boundary value of the third type enjoys a separate center discretization, and that the approaches to process discrete ADI scheme for i = 1 and i = M − 1 in the inner point vary widely. Therefore, only the discretization of the third type of boundary value and the processing mode of the discrete ADI format when i = 1 and i = M − 1 are discussed.
According to Equation (39), the total error R of the internal point is: whereT k 1,j is the virtual symmetric point of T k 1,j about x = 0, and is estimated to be [3]: According to Equations (62) and (63), we obtain: whereT k M−1,j is the virtual symmetric point ofT k M−1,j about x = X, and is estimated to be [3]: According to Equations (64)-(66), we have: whereT k i,M−1 is the virtual symmetric point of T k i,M−1 about r = R, and is estimated to be [3]: According to Equations (67)-(69), we have: For κ ∂T ∂r r=0 whereT k i,1 is the virtual symmetric point of T k i,1 about r = 0, and is estimated to be [3]: According to Equations (70) and (71), we have: If i = 1, the vector form of ADI with the backward Euler format is: −r +r 2 r j T k+1 The third kind of boundary discrete T k 0,j , T k M,j , T k i,M , T k i,0 , 0 ≤ i, j ≤ M are substituted into the above two equations to obtain the final discrete format for i = 1. The processing method for i = M − 1 performs the same as the one when i = 1. Therefore, the final matrix form of the difference scheme for the third kind of boundary value is: where, .
The other matrices, A, B, C, D, T, F, and Q are the same as the case of the first kind of boundary value.
The discretization of the third kind of boundary value is all the central difference quotient discretization, so there exists a normal number C 3 independent of h and τ, so that the total error R bz of boundary value discretization be: Therefore, Equations (61) and (78) represent that the normal number C 4 is independent of h and τ, so the total error R 1 of difference schemes for ADI with the backward Euler method of the third kind of boundary value is:

Difference Scheme for ADI with Crank-Nicolson Method
Based on the initial value condition, the boundary value condition and Equation (45), the difference scheme is obtained: The discrete processing of interior points and boundary values is the same as that of ADI with the backward Euler method for the third kind of boundary value. Therefore, the final matrix form of the difference scheme for the third kind of boundary value is: where, The other matrices, S, E, V, G, W, L, T, and ω are the same as the case of the first kind of boundary value.
Therefore, Equations (57) and (78) show that the normal number C 4 is independent of h and τ, so the total error R 1 of difference schemes for ADI with Crank-Nicolson method of the third kind of boundary value is: (87)

Stability of Difference Scheme Solutions
The following will prove that ADI format with the backward Euler method and ADI format with Crank-Nicolson method are unconditionally L 2 stable whether for the first kind of boundary value or the third kind of boundary value.

Definition 1. A function v(x) is defined on (−∞, +∞). If
where i = √ −1 is an imaginary unit and the above transformation is called the Fourier transform. Proof. From the total error R of the differential format for the first kind of boundary value and the total error R 1 of the differential format for the third kind of boundary value, the total error is O τ + h 2 .
Let λ = α·τ h 2 be the grid ratio and use the Fourier method to analyze the stability of Equation (25): Since the stability of homogeneous equations and non-homogeneous equations are consistent, the stability of homogeneous equations can be discussed. The difference equation of ADI with the backward Euler method corresponding to the homogeneous equation is: Substituting T k j,l = v k e ik 1 jh e ik 2 lh and e iθ = cos θ + i sin θ into the Equation (90), we get v k+1 =Ĝ(τ, k)v k , where the growth factorĜ(τ, k) is represented as: then, where k = (k 1 , k 2 ). Obviously, Ĝ (τ, k) ≤ 1 means that the difference scheme Equation (90) is unconditionally L 2 stable. Therefore, the difference Equation (25) is unconditionally L 2 stable. Proof. From the total error R (1) of the differential format for the first kind of the boundary value and the total error R (1) 1 of the differential format for the third kind of the boundary value, the total error is O τ 2 + h 2 .

Convergence of Difference Scheme Solutions
The following will prove that ADI format with the backward Euler method and ADI format with Crank-Nicolson method are convergent whether for the first kind of boundary value or the third kind of boundary value.

Definition 3.
For a sufficiently smooth function u, if the time step and the space step both approach 0, the truncation error of the difference equation approaches 0 for each node, it is said that the difference equation approximates the differential equation, that is, the difference equation is consistent with differential equations [2]. Theorem 3. If the difference equation satisfies the consistency condition and is stable according to the initial value, the difference solution converges to the solution of the original equation [2]. For the difference scheme Equation (25): T k i,j is expanded at node x i , r j , t k+1 for t, T k+1 i−1,j , T k+1 i+1,j and T k i−1,j , T k i+1,j are expanded at node x i , r j , t k+1 and x i , r j , t k for x, respectively,T k+1 i,j−1 , T k+1 i,j+1 and T k i,j−1 , T k i,j+1 are expanded at node x i , r j , t k+1 and x i , r j , t k for r, respectively. We have obtained: In the above equation, when h, τ → 0 , all the terms at the right-hand side of the above equation are close to 0, and the difference equation approaches to the original differential equation. Therefore, the difference Equation (25) T k i,j and T k+1 i,j is expanded at node x i , r j , t k+1 for t, T k+1 i−1,j , T k+1 i+1,j and T k i−1,j , T k i+1,j are expanded at node x i , r j , t k+1 and x i , r j , t k for x, respectively, T k+1 i,j−1 , T k+1 i,j+1 and T k i,j−1 , T k i,j+1 are expanded at node x i , r j , t k+1 and x i , r j , t k for r, respectively. We have obtained: In the above equation, when h, τ → 0 , all the terms at the right-hand side of the above equation are close to 0, and the difference equation approaches to the original differential equation. Therefore, the difference Equation (45) is consistent with the original equation. Proof. Order error is O τ 2 + h 2 , and according to theorems 3, 6, and 2, the difference scheme Equation (45) in the case of the first kind of boundary value and the third kind of boundary value is convergent. Using the difference scheme for ADI with the backward Euler method of the first kind of boundary value to solve the following definite solution problems, we can verify the convergence of the algorithm.

Example 1.
∂T ∂t T(x, r, 0) = e (x+r) , the exact solution of the above fixed solution problem is T(x, r, t) = e (x+r+t) .
and call x i , r j , t k as a network node.
The calculation formula of the maximum error E ∞ (h, τ) of the numerical solution is given by: The absolute error is the absolute value of the distinction between the exact solution and numerical solution. Table 1 shows the numerical solution, the exact solution, and the absolute error at some nodes when h = 1/32 and τ = 1/512. Table 1. Numerical solution, exact solution, and absolute error at some nodes (h = 1/32, τ = 1/512).

(x,r,t) Numerical Solution Exact Solution
Absolute Error Table 2 shows the maximum error of the numerical solution when the asynchronous step length is taken. It performs that when the space step is reduced to 1/2 and the time step is reduced to 1/4, the maximum error is approximately reduced by about 3/4.     Compare the longitudinal section and the cross section of Figures 2 and 3. If the similarity is high, it means that the numerical solution obtained by the algorithm is basically consistent with the numerical solution given by PDETOOL of MATLAB. It shows that the algorithm is more accurate.

Example of ADI with Crank-Nicolson Method
Using the difference scheme for ADI with Crank-Nicolson method of the first kind of boundary value to solve the following definite solution problems, we verify the convergence of the algorithm.
and call x i , r j , t k as a network node.
The calculation formula of the maximum error E ∞ (h, τ) of the numerical solution is given by: Table 3 shows the numerical solution, the exact solution, and the absolute error at some nodes when h = 1/32 and τ = 1/64.  Table 4 shows the maximum error of the numerical solution when the asynchronous step length is taken. It represents that when the space step is reduced to 1/2 and the time step is decreased to 1/2, the maximum error is approximately dropped to about 3/4.  Figure 6 shows the three-dimensional image of the approximate solution of the difference equation when h = 1/64 and τ = 1/128.   Compare the longitudinal section and the cross section of Figures 6 and 7. If the comparison is highly similar, the numerical solution obtained by the algorithm is basically consistent with the numerical solution given by PDETOOL of MATLAB, which shows that the algorithm is much accurate. Figure 8 shows the longitudinal section contrast between Figures 6 and 7.   From Figures 6-9, the numerical solution obtained by the ADI with Crank-Nicolson method for the first kind of boundary value is in good agreement with the numerical solution given by the PDETOOL of MATLAB in the same environment.
Corresponding to the comparison in the Tables 3 and 4 and Figures 6-9 that the ADI with Crank-Nicolson method has better convergence.

Laser Machining Simulation
We now use the difference scheme of the third kind of boundary value to solve the laser machining problem.
Example 3. Let a laser beam irradiate the surface of the material vertically under the ambient temperature of 293 K(20 • C). For simplicity, we consider that the input of laser energy acts on the material surface in the form of Gaussian heat flux. In this way, the center of the laser can reach the heat peak, and the temperature of the material boundary tends to the initial temperature. Therefore, the laser radius is equal to the material radius. Table 5 shows the thermal property parameters of 316 stainless steel. 2.
In the laser irradiation area Γ, the moving laser beam is loaded through the boundary conditions of the surface heat source: 3.
The boundary outside the laser irradiation area Γ is in contact with air, and the boundary conditions are as following: 4.
The initial temperature of the substrate is the ambient temperature, that is, the initial conditions of the substrate are: Cylindrical coordinates should be adopted to this three-dimensional problem. Since the temperature distribution along the depth is symmetric after the laser beam irradiates the surface of the material, it can be converted into a two-dimensional problem with coordinates (r, x). The heat conduction format in cylindrical coordinates is rewritten as: where α = κ ρc is the thermal diffusivity. Using ADI with the backward Euler method and ADI with Crank-Nicolson method to solve the problem on MATLAB, we can obtain Figures 10 and 11.    The laser source is a function of r, and the temperature is symmetrically distributed along the radius of the material. Therefore, the temperature becomes lower in the x direction, distributed along the laser source function in the r direction. Table 6 is the maximum error and the maximum relative error between the numerical solutions of the ADI with the backward Euler method and the exact solutions of ADI with Crank-Nicolson method when the asynchronous step length is taken. The maximum relative error is the maximum error divided by the corresponding exact solution.
From Table 6, when τ/h gets smaller, the maximum error and the maximum relative error of the two algorithms gets smaller. The smaller τ/h 2 , the less the maximum error and the maximum relative error of the two algorithms when τ/h comes the same.
Use the PDETOOL of MATLAB to simulate Example 3. The results are shown in Figure 12. The approximate solution in Figures 10 and 11 is compared with the cross and the longitudinal sections of the numerical solution in Figure 12, as shown in Figures 13 and 14.    Figure 14 shows the cross-sectional contrast among Figures 10-12. From Figures 10-14, both numerical solutions obtained by the ADI with the backward Euler method and with Crank-Nicolson method for the third kind of boundary value are in good agreement with the numerical solutions given by the PDETOOL of MATLAB in the same environment within a certain small error, which shows that these methods are feasible and effective in laser processing.

Example 4.
To be compared with Example 3, we take the laser radius as 0.0008 m so that the laser radius is smaller than the material radius under the same environmental conditions.
(2) The input of laser energy acts on the material surface in the form of Gaussian heat flux, specifically: Other parameters and environmental conditions behavior the same as Example 3.
Using the ADI with Crank-Nicolson method to solve the problem on MATLAB, we can obtain Figure 15.  Use the PDETOOL of MATLAB to simulate Example 4. The results are shown in Figure 16.
From the results of Figures 15 and 16, it shows the correctness of the algorithm. We consider combining the βq(r) of Example 4 and the rest of the environmental conditions of Example 3 to form a new problem. This new problem and Example 4 vary from the laser irradiation area. Using ADI with Crank-Nicolson method to solve this new problem on MATLAB, we can attain Figure 17.   From Figures 15 and 17, the maximum error and the maximum relative error obviously appear on the boundary. Since the maximum error and the maximum relative error respectively is 8.5667 K and 0.0284, and the melting point of the material is relatively high, the temperature distinction in the interior point is negligible in Figures 15 and 17.

Conclusions
In this article, the effective computation of the heat distribution of the material is studied when the laser beam is irradiated on the section of the cylinder material, where the laser beam is expressed as a Gaussian distribution. The mathematical model-threedimensional heat conduction equation, is converted into a two-dimensional parabolic equation. Based on the symmetry of the heat distribution, the three-dimensional equation in the rectangular coordinate system can be changed to the simplified two-dimensional equation in the cylindrical coordinate system, which raises the work efficiency, but also saves the storage space.
Subsequently, to solve the simplified equation, the unconditionally stable ADI scheme is developed with the classic backward Euler method or Crank-Nicolson method in the cylindrical coordinate system, and then the simulation results are obtained, where the first kind of boundary value condition or the third kind of boundary value condition is considered. Further, these two methods have been proved to be unconditionally L 2 stable and convergent, and their accuracies are verified by Examples 1 and 2. Examples 3 and 4 demonstrate that these two methods are feasible, stable, and efficient to solve the laser processing problem model. Example 4 also verifies that under certain conditions and errors, the entire laser irradiation can be used to replace the partial laser irradiation, thereby simplifying the boundary heat dissipation problem.
Comparison between the results obtained by these two methods and the results obtained by PDETOOL of MATLAB in each example, shows that our results are more convincing, and our treatment has two advantages. Firstly, self-programming can use different methods to discretize distinct problems, which improves the flexibility of the algorithm, saves running space and time and raises efficiency. Secondly, in the absence of an exact solution, our approach can check whether the algorithm is correct or not by comparing with the numerical solution of PDETOOL.
The three-dimensional heat conduction problem model of the article is also universal. The material irradiated by the laser is not necessarily cylindrical. As long as a circle is set at the center of the laser spot, the problem model in this article is still applicable. The only difference is the processing of the boundary, that is, the expression and discretization of the boundary conditions may need to be processed separately.