Abstract
This paper focuses on efficiently numerical investigation of two-dimensional heat conduction problems of material subjected to multiple moving Gaussian point heat sources. All heat sources are imposed on the inside of material and assumed to move along some specified straight lines or curves with time-dependent velocities. A simple but efficient moving mesh method, which continuously adjusts the two-dimensional mesh dimension by dimension upon the one-dimensional moving mesh partial differential equation with an appropriate monitor function of the temperature field, has been developed. The physical model problem is then solved on this adaptive moving mesh. Numerical experiments are presented to exhibit the capability of the proposed moving mesh algorithm to efficiently and accurately simulate the moving heat source problems. The transient heat conduction phenomena due to various parameters of the moving heat sources, including the number of heat sources and the types of motion, are well simulated and investigated.
1. Introduction
Heat conduction phenomena of material involving moving heat sources, which have attracted increasing attention by scientists and engineers in the past few decades, have been studied in a wide range of fields, such as welding, cutting, drilling, laser hardening/forming, plasma spraying, heat treating of metals, manufacturing of electronic components, and even firing a gun barrel, solid propellant burning, and dental treatment, see e.g., [1–5] and references therein. The most important physical quantity of interest for such practical applications is the temperature field of the medium, which is usually modeled by the heat conduction equation with time-dependent localized source terms for moving heat sources. Once the temperature field is obtained, many other thermophysical properties of material, including metallurgical microstructures, thermal stress, residual stress, and part distortion, could be subsequently determined [6–10]. It is therefore particularly important to precisely and efficiently predict the dynamic variation of the temperature field around the moving heat sources during these engineering processes.
In order to investigate the temperature field and the related thermal properties of the problem with moving heat sources, numerous methods, in either analytical or numerical approach, have been developed, since the 1930s, when the pioneering work of Rosenthal was proposed for the analytical solution of a simplified moving heat source problem [11]. Although analytical methods are still popular nowadays [12], they are usually only available for simple situations such as the quasistationary problem of a single heat source moving along a straight line with a constant speed. In comparison to analytical methods, numerical methods could only provide results approximately within an acceptable error tolerance, but they are more flexible to deal with the complicated yet practical situations such as the transient problem of multiple heat sources moving in a complex geometry of the material with time-dependent speeds [3]. However, most of numerical studies, regardless using meshless methods [13, 14] or mesh-based methods such as the finite element method [6, 10], were concerned about problems involving only a heat source moving along a straight line with a constant speed, or multiple heat sources moving along parallel straight lines with the same constant speed. Apart from these, the technique of moving coordinate system, such that the heat source is stationary in the new coordinate system, is often introduced in both analytical and numerical analyses of the quasistationary problem [1, 13]. Nevertheless, it is obvious that this technique is limited and not applicable for problems subjected to multiple moving heat sources with different velocities and trajectories.
It is well known that the moving heat source might be imposed on the surface or inside of material [2], which follows that the resulting mathematical model would contain a source term in the boundary condition or the governing heat conduction equation, respectively. Depending on the practical applications, the moving heat source can be modeled as a point, line, or plane source with various geometries, such as square, circle, semiellipsoidal and double ellipsoidal [1, 15, 16]. No matter what kind of the moving heat source, its energy is always highly concentrated in a time-dependent localized domain. It turns out that the resulting temperature of material would change drastically in the localized region around the moving heat source. Consequently, it is obvious that a significant improvement in efficiency could be achieved, if an adaptive mesh method, which concentrates a number of mesh points dynamically in the local regions of rapid variation of the temperature, is employed to solve the problem with the same accuracy as the fixed mesh method.
The moving mesh method [17, 18] is one of the popular adaptive methods and has been successfully applied to various problems that contain time-dependent localized singularities [19–21]. It usually tries to find a time-dependent one-to-one coordinate transformation between the physical domain and the computational domain, by solving an additional system of moving mesh partial differential equation (MMPDE), which equidistributes a certain monitor function of the physical solution [22, 23]. The original physical equation would be subsequently transformed into the computational domain and then be solved by the standard uniform mesh method. For more details of the moving mesh method, one is referred to [17, 23, 24]. Up to now, the moving mesh method has been exhibited to work well for problems with moving heat source in a one-dimensional (1D) case [25, 26]. Yet the application of the moving mesh method to the problem of moving heat source in multidimensional case is still immature.
Based on the above observations, this paper is concerned about the efficiently numerical study of two-dimensional (2D) heat conduction problems involving multiple moving heat sources by the moving mesh method. The Gaussian point heat source, that is imposed on the inside of material and allowed to move along any specified curve with a time-dependent velocity, is taken for all heat sources as an example of the model problem. A simple moving mesh method, which generates the 2D moving mesh dimension by dimension from 1D MMPDE with an appropriately defined monitor function, is developed. The transient heat conduction phenomena due to various parameters of the moving heat sources, such as the number of heat sources and the types of motion, are then investigated with the proposed moving mesh method. Since only two additional 1D systems are required to be solved, the resulting moving mesh method is easy to be implemented and turns out to be very efficient to give satisfactory results.
The rest of the paper is outlined as follows. In Section 2, the mathematical model of the 2D heat conduction problem with multiple moving heat sources is briefly introduced. The detailed formulation of the moving mesh method for the model problem is described in Section 3. Numerical experiments are presented to show the efficiency of the proposed moving mesh method in Section 4, where heat conduction phenomena are also investigated in detail. Finally, some conclusions are given in the last section.
2. Mathematical Model
In a thin rectangular plate made of homogeneous material, heat flow can be simplified to be viewed as a two-dimensional flow. Let the plate occupy domain , where and are the length and width of the plate, respectively. Suppose the plate is initially at room temperature denoted by and is heated by several moving heat sources at time , as shown in Figure 1. Then using to represent the temperature at position and time , the evolution of the temperature in the plate can be described by the following two-dimensional heat conduction equation: where , , and are the material density, the heat capacity, and the thermal conductivity, respectively. In the current investigation, these quantities are assumed to be constant independent of the position and temperature. The right-hand side of (1) represents the heat source term, where is the number of heat sources and is the volumetric heat generation rate of the th heat source.
Depending on the physical nature of the problem, a moving heat source can be roughly classified into three types, namely, the point, line, and plane heat source. All of them concentrate high power in a time-dependent localized region and can be well modeled by a Dirac delta function [1, 2, 8, 12]. However, the singularity of delta function introduces additional difficulties especially for numerical simulation of practical engineering applications. Consequently, a well-defined smooth function such as the localized Gaussian distribution function is usually introduced to replace the delta function when researchers study the problem from numerical approaches [6, 10, 13, 14]. In this paper, we are mainly interested in the heat conduction due to moving Gaussian point heat sources, which takes the form for the th heat source. Here, is the effective heating radius of the th heat source, and is the maximum heat flux at the center of the corresponding heat source, whose moving trajectory is given by .
To complete the description of the problem, it remains to give the initial condition at time and the boundary condition throughout the simulation time . Obviously, we have the initial condition from the previous assumption. For the boundary condition, it is convenient to divide the boundary of the plate into two parts, i.e., , and let where and are the prescribed temperature and heat flux, respectively, and is the unit outward normal vector. In other words, the Dirichlet boundary condition is applied on , while the Neumann boundary condition is applied on .
At last, it is noted that the above 2D model is also able to describe the temperature evolution with the moving line heat source, as shown in [1, 2].
3. Formulation of the Numerical Method
This section is devoted to illustrate the details of the moving mesh method to solve the model problem (1)–(3). We first give a brief review of the 1D moving mesh partial differential equation. Based on it, a strategy of 2D moving mesh generation is introduced. The discretization of the model equations on the resulting moving mesh, together with the final algorithm of numerical simulation, will then be presented.
3.1. 1D Moving Mesh Partial Differential Equation
Let and denote the physical and computational coordinates, respectively. A time-dependent one-to-one coordinate transformation between the physical domain and the computational domain, which are without loss of generality assumed to be and , respectively, is denoted by with and . For a uniform mesh on the computational domain, given by with , a time-dependent mesh on the physical domain can be correspondingly obtained by setting for all . Therefore, in order to find an adaptive physical mesh that dynamically concentrates mesh points in regions of interest, e.g., the regions of a rapid variation of the solution, it is equivalent to find a suitable coordinate transformation according to some special measure of the solution.
Based on the equidistribution principle, such a transformation can be obtained by solving the following equation [17, 22]: with boundary conditions and . Here, is a user-defined function of the solution to control the concentration of the mesh. It is called the monitor function or the mesh density function in the theory of the moving mesh method and will be given specifically in Section 3.4 for our numerical experiments.
In practice, the quasistatic equation (5) is usually relaxed by adding terms involving the mesh speed . The resulting equation is referred to as a moving mesh partial differential equation (MMPDE). Among the various MMPDEs proposed over the past few decades, we would like to utilize the so-called MMPDE6 [22] in the present work, since it has been shown to work well for the moving heat source problem [25, 26]. The MMPDE6 reads where is a positive parameter for adjusting the response time of mesh movement to the change of the monitor function . With boundary conditions and , the adaptive physical mesh would be updated at the moment by solving the linear system derived from the finite difference discretization of MMPDE6, that is, for , where is the time step length, is the numerical approximation of the th mesh point at time , and with is the discrete monitor function on the th mesh point at time . Nevertheless, it is pointed out that the MMPDE6 could also be solved by the MATLAB package called MMPDElab [23].
3.2. 2D Moving Mesh Generation
A complete two-dimensional MMPDE and the resulting moving mesh method, as can be seen in [17], are in some sense a little complicated and not easy to use. On the other hand, an adaptive rectangular mesh on the physical domain generated by 1D mesh strategy is obviously much simpler and has also been successfully applied to reaction-diffusion equations of quenching type, see e.g. [19, 27]. Accordingly, we shall follow the later approach to generate the adaptive rectangular mesh on the physical domain via 1D MMPDE dimension by dimension in this paper.
To be specific, let the time-dependent one-to-one coordinate transformation between 1D domains and be still denoted by with and . Given a uniform mesh on the domain with for , a time-dependent mesh on the domain could be obtained by setting for all . Similarly, by introducing a time-dependent one-to-one coordinate transformation with and between 1D domains and , a time-dependent mesh on the domain could be obtained by , where with is the uniform mesh on the domain . Then a time-dependent rectangular mesh on the physical domain would be generated by setting the mesh point to be for all and .
As stated in the previous subsection, both and can be determined from 1D MMPDE6 by utilizing appropriate monitor functions and , respectively, where and are functions of the 2D solution , and will be specified in Section 3.4.
Obviously, the above strategy of 2D moving mesh generation is very efficient and easy to be implemented, since only two one-dimensional linear systems are need to be solved.
3.3. Discretization on the Moving Mesh
It is now ready to introduce the discretization of the model equations (1)–(3) on the 2D rectangular moving mesh using the central finite difference method.
Using the time-dependent coordinate transformations and between the physical coordinates , and the computational coordinates , , any function of , , and can be expressed as a function in terms of , , and , that is, . By the chain rule, it follows that
In order to distinguish the two partial derivatives with respect to in the above expression, the notation , similar to the notation of mesh speed , is introduced for the first one, i.e., , and the other one is simplified to the original notation without causing confusion. Then, in the computational coordinates , and , the original physical equation (1) becomes where is the thermal diffusivity and .
The above equation can be discretized using the second-order central finite difference method on the uniform computational mesh with and . This subsequently yields a system of ordinary differential equations of the form where , , , and
Using the Crank-Nicolson method for temporal discretization, a full discretization, which has second-order time accuracy, can be obtained by for and . In the above equation, is the numerical approximation of the temperature at at time , equivalently at of the physical domain at time . As for , , and , they are numerical approximations of , , and , respectively, and computed by substituting all time-dependent quantities with the corresponding numerical approximations in (11). Similarly, , , , and are corresponding numerical approximations at time .
Supplemented with appropriate discretization of boundary condition (3), the linear system (12) can then be solved for all . Let us take the left boundary where or equivalently as an example. If on the left boundary it is subjected to the Dirichlet boundary condition, we shall directly set for all . Alternatively, if on the left boundary it is subjected to the Neumann boundary condition, which reduces to we shall take the following discretization: for all , to make sure the discretization of the boundary condition is also second-order accuracy.
3.4. Final Algorithm and the Monitor Function
Now, we are in a position to describe the whole numerical algorithm that simulates the moving heat source problem with the moving mesh method. It is evident that the full discretization, including the system of the discretization (12) and the discretization of two 1D MMPDE6 for and , respectively, is coupled together via the monitor functions and the physical mesh. A simple decouple strategy is adopted in the present algorithm, that is, the mesh equation and the physical equation are solved alternately one by one. A flowchart of the final moving mesh algorithm is then presented in Algorithm 1. However, to close this section, it remains to give the details of the monitor functions and .
It is well known that the monitor function plays an important role to the success of the moving mesh method [17]. One of the popular choices is the arc-length monitor function, which is aimed at equidistributing the arc-length of the solution curve between each two adjacent mesh points. As a result, it usually works well and is able to concentrate the mesh points in the local regions of a large derivative of the solution. Additionally, if there are local regions with large curvature of the solution, then the curvature monitor function might be a good candidate.
For the moving heat source problem, it is easy to show that there are not only local regions with large derivatives of the solution but also local regions, e.g., the neighborhood of the point heat source, where the curvature of the solution is large and the derivative is close to . In view of these, a linear combination of the arc-length monitor function and the curvature monitor function, which reads is employed in our numerical experiments. Here, is a 1D function defined later by a certain average of the 2D temperature with respect to , and is the weight of the arc-length monitor function. Applying the central finite difference method to (16), one can obtain on by where . Apparently, it is enough to give in the computation of . Taking the whole 2D temperature field into consideration, the value of may be defined by
Furthermore, it is pointed out in [17] that the smoothness of the monitor function may affect the stability and quality of the moving mesh. Consequently, is smoothed in our simulations by the strategy proposed in [28], i.e., where and are two smoothing parameters, given by and currently. Following the same approach, we can get by replacing , , and in the right-hand side of (17), respectively, with , , and Then, would be smoothed with the same strategy of (19).
| ||||||||||||||||||||||
4. Numerical Experiments and Discussion
Several numerical experiments are carried out in this section to show the capability to efficiently and accurately simulate moving heat source problems with the proposed algorithm, which is implemented in MATLAB (Release 2016a, The MathWorks, Inc., Natick, Massachusetts, MA, USA). Heat conduction phenomena in the plate due to the number of moving point heat sources, the types of motion, and some other properties are also investigated in detail.
Throughout the simulation, the units presented in Table 1 are employed for the involved physical variables, and the plate is assumed to be homogeneous with the material density , the heat capacity , and thermal conductivity . The room temperature is adopted for both the initial temperature and the boundary temperature . When the Neumann boundary condition is taken into account, the boundary heat flux would be . Moreover, the time step length is given by , the parameter in MMPDE6 takes the value of , and the weight in the monitor function is set to be , if they are not explicitly pointed out below. For the rest of the parameters, they will be specified for each experiment individually.
4.1. A Heat Source Moving along a Straight Line
The first experiment focuses on the case that the plate is subjected to a single Gaussian point heat source, which moves along the -axis with a constant speed. A lot of research, including both numerical and analytical studies, can be found in the literature for this case. Here, the same problem settings as in [13, 14] are considered. To be specific, the plate has the length of and the width of . The Dirichlet boundary condition is applied on the left boundary of the plate, while the rest of the boundaries of the plate are subjected to the Neumann boundary condition. The moving Gaussian point heat source, defined by the effective radius of and the maximum heat flux of , is initially at the center of the right boundary and moves from right to left along -axis with a constant speed of . It follows that and .
The test is simulated on the moving mesh with several different values of and . Numerical results are subsequently compared with the solutions obtained from the discretization (12) on the uniform mesh with various and . As presented in Figure 2, for the temperature profile along the heat source moving path, i.e., -axis with , at , it can be seen that the solution on the moving mesh with and is much more accurate than the solution on the uniform mesh with the same and . In fact, it is comparable to the results in the uniform mesh with and . In each time step, a single linear system of order is required to be solved for the uniform mesh algorithm, whereas for the moving mesh algorithm, three linear systems of order , , and , respectively, are required to be solved. It follows obviously that the proposed moving mesh algorithm is able to give the solution within the same accuracy more efficiently than the related uniform mesh algorithm.
(a)
(b)
(c)
The transient 2D temperature field obtained by the moving mesh algorithm with and , together with the corresponding physical mesh, is presented in Figure 3, at time instances (a), 25, 35, and 45 (d), respectively. It turns out that these temperature fields show a good agreement with the results reported in [13, 14]. It is also found that during the simulation, the physical mesh could be adjusted successfully and dynamically according to the temperature field, so that the algorithm always concentrates a number of mesh points in regions of interest as the monitor function indicated.
(a)
(b)
(c)
(d)
At last, it is worth mentioning that the peak temperature occurs near the rear of the moving heat source, rather than the exact position of the heat source, as can be observed in Figure 2. This is not surprising and can be understood by noting that the moving heat source is always exposed to a much cooler position, and the temperature near the rear of the heat source may continue to increase if the heat does not spread out in time. In addition, similar phenomena have been observed from the results reported in [14].
4.2. A Heat Source Moving along a Circle
The second experiment considers the case that a square plate with side length is subjected to a single Gaussian point heat source, which moves along a circle of radius with a constant speed in a counterclockwise direction. Specifically, the heat source has the effective radius of and the maximum heat flux of . Its moving path is set to be and . Additionally, all boundaries of the plate are assumed to satisfy the Dirichlet boundary condition.
This experiment is simulated by the moving mesh algorithm with and the weight in the monitor function to be . The transient 2D temperature field at time instances (a), 2, 3, and 4 (d), respectively, as well as the corresponding physical mesh, is depicted in Figure 4, where the pink circle represents the heat source moving path. As can be seen from Figure 4, the moving mesh algorithm still successfully concentrates enough mesh points in regions of interest as the monitor function indicated.
(a)
(b)
(c)
(d)
As a result, the proposed algorithm is also able to be employed to investigate heat conduction phenomena for this case accurately with a small number of and . Thus, a great improvement in efficiency can be obtained by the proposed algorithm.
After a long time simulation, a quasistationary temperature field can be achieved. As shown in Figure 5, it would be stationary in the moving coordinate system that attaches to the moving heat source. Similar results can also be found in [29].
4.3. Multiple Heat Sources Moving along Straight Lines
Now let us go to investigate the heat conduction phenomena of the plate subjected to multiple moving heat sources. Three cases, that is, two heat sources moving along -axis in opposite directions, two heat sources moving along two intersecting straight lines, and three heat sources moving along three straight lines parallel to -axis, are considered below. In all cases, the Dirichlet boundary condition is adopted for the left boundary of the plate, while the Neumann boundary condition is adopted for the rest of the boundaries of the plate. The all involved heat sources are assumed to be Gaussian point heat source with the effective radius to be and the maximum heat flux to be , except for the last case where .
For the first case, the size of the plate is set to be and . The two heat sources are suddenly imposed on the position , respectively, at the initial time, and then move along -axis in opposite directions with a constant speed of . The resulting moving paths are and .
Obviously, the two heat sources will meet each other at time instance . The simulation is performed by the proposed moving mesh algorithm with and . The corresponding transient 2D temperature field as well as the physical mesh are presented in Figure 6, for time instances (a), , , and (d), respectively. Additionally, the 1D temperature profiles along the heat source moving path at time instances , , , , , and are given in Figure 7. It can be seen that the physical mesh moves adaptively according to the monitor function of the temperature field, in which there exists a peak following each heat source. As the two heat sources approach each other, two peaks would merge into a single peak, causing the peak temperature to increase rapidly to a high level near . Then two peaks are separated as the two heat sources move away from each other, and the peak temperature subsequently decreases to the normal level around .
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(e)
(f)
For the second case, the plate is square with side length . The heat source moving paths are set to be . That is, the two heat sources are initially at the position , and move along the straight lines , respectively, with the constant speed of . Thus, they will meet each other at the original point at time instance . The transient 2D temperature field and the corresponding physical mesh, obtained by the moving mesh algorithm with , are plotted in Figure 8 for time instances (a), , , and (d), respectively. Similar phenomena could be observed as the previous case.
(a)
(b)
(c)
(d)
For the last case, the plate is the same as the first case, i.e., and . The moving paths of the three heat sources are set to be , , and , which follows that the three heat sources are initially at the right boundary and move from right to left along horizontal lines with the same constant speed of . The resulting transient 2D temperature field and the corresponding physical mesh by the moving mesh algorithm with and are shown in Figure 9 for time instances (a), , , , and (e), respectively. Again, the physical mesh could be adjusted successfully according to the monitor function of the temperature field. Since the heat sources move much faster than other cases, the peak temperature in this case is smaller than the one in the previous cases.
(a)
(b)
(c)
(d)
(e)
5. Conclusions
A simple moving mesh algorithm has been developed to numerically solve the 2D model equations of moving heat source problems with Gaussian point heat sources. In the present algorithm, only two additional 1D mesh equations are required to be solved for each time step. However, it is found that the physical mesh could successfully and dynamically concentrate a number of mesh points in regions of interest as the monitor function indicated. Therefore, the proposed algorithm is able to simulate the moving heat source problem very accurately and efficiently. Heat conduction phenomena of the rectangular plate subjected to moving Gaussian point heat sources with various types of motion, including moving along straight lines and a circular trajectory, have then been numerically investigated. Numerical results validate the accuracy and efficiency of the proposed algorithm, which shows that the proposed moving mesh algorithm is a promising approach for such moving heat source problems.
Finally, the extension of the proposed moving mesh algorithm to other localized heat source models, such as Dirac delta point heat source and plane heat source, is ongoing and would be presented elsewhere soon. The full 3D simulation of the moving heat source problem with the moving mesh method will also be studied in the future work.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The authors would like to express their sincere appreciation to the editor and the anonymous referees for their valuable comments and suggestions. This work was supported by the National Natural Science Foundation of China (grant number 11601229) and the Jiangsu Province Natural Science Foundation of China (grant number BK20160784).