# Heat Conduction Simulation of 2D Moving Heat Source Problems Using a Moving Mesh Method

Authors: Zhicheng Hu; Zhihui Liu
Journal: Advances in Mathematical Physics
 (2020)

---

## 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.


---

## Body

## 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 Section2, 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Ω=x,y:−Lx/2≤x≤Lx/2,−Ly/2≤y≤Ly/2, where Lx and Ly are the length and width of the plate, respectively. Suppose the plate is initially at room temperature denoted by T0 and is heated by several moving heat sources at time t>0, as shown in Figure 1. Then using Tx,y,t to represent the temperature at position x,y and time t, the evolution of the temperature in the plate can be described by the following two-dimensional heat conduction equation:
(1)ρc∂T∂t−k∂2T∂x2+∂2T∂y2=∑l=1qglx,y,t,x,y∈Ω,t>0,where ρ, c, and k 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 q is the number of heat sources and glx,y,t is the volumetric heat generation rate of the lth heat source.Figure 1
Sketch of the rectangular plate with a moving point 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
(2)glx,y,t=Qlexp−x−αlt2+y−βlt2rl2,for the lth heat source. Here, rl is the effective heating radius of the lth heat source, and Ql is the maximum heat flux at the center of the corresponding heat source, whose moving trajectory is given by αlt,βlt.To complete the description of the problem, it remains to give the initial condition at timet=0 and the boundary condition throughout the simulation time t. Obviously, we have the initial condition Tx,y,0=T0 from the previous assumption. For the boundary condition, it is convenient to divide the boundary of the plate into two parts, i.e., ∂Ω=Γ1∪Γ2, and let
(3)T=T¯onΓ1,n·k∇T=q¯onΓ2,where T¯ and q¯ are the prescribed temperature and heat flux, respectively, and n is the unit outward normal vector. In other words, the Dirichlet boundary condition is applied on Γ1, while the Neumann boundary condition is applied on Γ2.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

Letx 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 a,b and 0,1, respectively, is denoted by
(4)x=xξ,t,ξ∈0,1,with x0,t=a and x1,t=b. For a uniform mesh on the computational domain, given by ξj=j/N with j=0,1,⋯,N, a time-dependent mesh on the physical domain can be correspondingly obtained by setting xjt=xξj,t for all j. 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 xξ,t 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]:
(5)∂∂ξMxξ,t,t∂x∂ξ=0,with boundary conditions x0,t=a and x1,t=b. Here, Mx,t 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 ẋξ,t=∂/∂txξ,t. 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
(6)∂2ẋ∂ξ2=−1τ∂∂ξM∂x∂ξ,where τ is a positive parameter for adjusting the response time of mesh movement to the change of the monitor function Mx,t. With boundary conditions x0t=a and xNt=b, 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,
(7)xj+1n+1−2xjn+1+xj−1n+1−xj+1n−2xjn+xj−1nΔtn=−1τMj+1/2nxj+1n+1−xjn+1−Mj−1/2nxjn+1−xj−1n+1for j=1,2,⋯,N−1, where Δtn=tn+1−tn is the time step length, xjn≈xjtn is the numerical approximation of the jth mesh point at time tn, and Mj+1/2n=Mj+1n+Mjn/2 with Mjn=Mxjn,tn is the discrete monitor function on the jth mesh point at time tn. 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−Lx/2,Lx/2 and 0,1 be still denoted by x=xξ,t with x0,t=−Lx/2 and x1,t=Lx/2. Given a uniform mesh on the domain 0,1 with ξi=i/Nx for i=0,1,⋯,Nx, a time-dependent mesh on the domain −Lx/2,Lx/2 could be obtained by setting xit=xξi,t for all i. Similarly, by introducing a time-dependent one-to-one coordinate transformation y=yη,t with y0,t=−Ly/2 and y1,t=Ly/2 between 1D domains −Ly/2,Ly/2 and 0,1, a time-dependent mesh on the domain −Ly/2,Ly/2 could be obtained by yjt=yηj,t, where ηj=j/Ny with j=0,1,⋯,Ny is the uniform mesh on the domain 0,1. Then a time-dependent rectangular mesh on the physical domain Ω would be generated by setting the mesh point to be xit,yjt for all i and j.As stated in the previous subsection, bothxit and yjt can be determined from 1D MMPDE6 by utilizing appropriate monitor functions Mx,t and Gy,t, respectively, where Mx,t and Gy,t are functions of the 2D solution Tx,y,t, 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 transformationsx=xξ,t and y=yη,t between the physical coordinates x, y and the computational coordinates ξ, η, any function of x, y, and t can be expressed as a function in terms of ξ, η, and t, that is, fx,y,t=fxξ,t,yη,t,t. By the chain rule, it follows that
(8)∂f∂ξ=∂f∂x∂x∂ξ,∂f∂η=∂f∂y∂y∂η,∂f∂tξ,ηfixed=∂f∂tx,yfixed+∂f∂xẋξ,t+∂f∂yẏη,t.In order to distinguish the two partial derivatives with respect tot in the above expression, the notation ḟ, similar to the notation of mesh speed ẋ, is introduced for the first one, i.e., ḟ=∂f/∂tξ,ηfixed, and the other one is simplified to the original notation ∂f/∂t without causing confusion. Then, in the computational coordinates ξ, η∈0,1 and t>0, the original physical equation (1) becomes
(9)Ṫ−∂T∂ξ/∂x∂ξẋ−∂T∂η/∂y∂ηẏ−μ∂∂ξ∂T∂ξ/∂x∂ξ/∂x∂ξ+∂∂η∂T∂η/∂y∂η/∂y∂η=g~ξ,η,t,where μ=k/ρc is the thermal diffusivity and g~ξ,η,t=1/ρc∑l=1qglxξ,t,yη,t,t.The above equation can be discretized using the second-order central finite difference method on the uniform computational meshξi,ηj with i=0,1,⋯,Nx and j=0,1,⋯,Ny. This subsequently yields a system of ordinary differential equations of the form
(10)Ṫi,j−Ai,jtẋi−Bi,jtẏj−Li,jt=0,i=1,2,⋯,Nx−1,j=1,2,⋯,Ny−1,where Ti,j=Ti,jt=Tξi,ηj,t, xi=xit, yj=yjt, and
(11)Ai,jt=Ti+1,jt−Ti−1,jtxi+1t−xi−1t,Bi,jt=Ti,j+1t−Ti,j−1tyj+1t−yj−1t,Li,jt=μ2xi+1t−xi−1tTi+1,jt−Ti,jtxi+1t−xit−Ti,jt−Ti−1,jtxit−xi−1t+2yj+1t−yj−1tTi,j+1t−Ti,jtyj+1t−yjt−Ti,jt−Ti,j−1tyjt−yj−1t+g~ξi,ηj,t.Using the Crank-Nicolson method for temporal discretization, a full discretization, which has second-order time accuracy, can be obtained by(12)Ti,jn+1−Ti,jnΔtn−Ai,jn+Ai,jn+12xin+1−xinΔtn−Bi,jn+Bi,jn+12yjn+1−yjnΔtn−Li,jn+Li,jn+12=0,for i=1,2,⋯,Nx−1 and j=1,2,⋯,Ny−1. In the above equation, Ti,jn≈Ti,jtn is the numerical approximation of the temperature at ξi,ηj at time tn, equivalently at xin,yjn of the physical domain at time tn. As for Ai,jn, Bi,jn, and Li,jn, they are numerical approximations of Ai,jtn, Bi,jtn, and Li,jtn, respectively, and computed by substituting all time-dependent quantities with the corresponding numerical approximations in (11). Similarly, Ti,jn+1, Ai,jn+1, Bi,jn+1, and Li,jn+1 are corresponding numerical approximations at time tn+1.Supplemented with appropriate discretization of boundary condition (3), the linear system (12) can then be solved for all Ti,jn+1. Let us take the left boundary where x=−Lx/2 or equivalently ξ=0 as an example. If on the left boundary it is subjected to the Dirichlet boundary condition, we shall directly set
(13)T0,jn+1=T¯,for all j. Alternatively, if on the left boundary it is subjected to the Neumann boundary condition, which reduces to
(14)−k∂T∂x=q¯,we shall take the following discretization:
(15)T1,jn+1−T0,jn+1x1n+1−x0n+1+T2,jn+1−T0,jn+1x2n+1−x0n+1−T2,jn+1−T1,jn+1x2n+1−x1n+1=−q¯kfor all j, 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 xin+1 and yjn+1, 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 Min and Gjn.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 to0. In view of these, a linear combination of the arc-length monitor function and the curvature monitor function, which reads
(16)Mx,t=θ1+∂u∂x2+1−θ1+∂2u∂x224,is employed in our numerical experiments. Here, u=ux,t is a 1D function defined later by a certain average of the 2D temperature Tx,y,t with respect to y, and θ is the weight of the arc-length monitor function. Applying the central finite difference method to (16), one can obtain Min on xin by
(17)Min=θ1+ui+1n−ui−1nxi+1n−xi−1n2+1−θ1+2xi+1n−xi−1nui+1n−uinxi+1n−xin−uin−ui−1nxin−xi−1n24,where uin=uxin,tn. Apparently, it is enough to give uin in the computation of Min. Taking the whole 2D temperature field into consideration, the value of uin may be defined by
(18)uin=1Ny+1∑j=0NyTi,jn.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, Min is smoothed in our simulations by the strategy proposed in [28], i.e.,
(19)Min≔∑k=i−νi+νMkn2γ1+γ∣k−i∣/∑k=i−νi+νγ1+γ∣k−i∣,where γ>0 and ν≥0 are two smoothing parameters, given by γ=2 and ν=2 currently. Following the same approach, we can get Gjn by replacing i, xin, and uin in the right-hand side of (17), respectively, with j, yjn, and
(20)u~jn=1Nx+1∑i=0NxTi,jn.Then, Gjn would be smoothed with the same strategy of (19).Algorithm 1. Flowchart of the moving mesh algorithm for moving heat source problem.
Input: The end time tend, initial physical mesh (xi0, yj0) and initial temperature field Ti,j0.
Output: The final physical mesh (xin, yjn) and the corresponding temperature field Ti,jn.
1 Letn=0 and tn=0;
2while tn<tend do
3  Determine the time step∆tn;
4  Compute the 1D monitor functionsMin on xin for all i, and Gjn on yjn for all j, based on the current physical mesh (xin, yjn) and the corresponding temperature field Ti,jn;
5  Solve two linear systems of discretization of 1D MMPDE6 withMin and Gjn, respectively, to get two new 1D mesh xin+1 and yjn+1;
6  Construct the new physical mesh (xin+1, yjn+1);
7  Solve the system of discretization (12) to get the new temperature field Ti,jn+1;
8  Lettn+1=tn+∆tn and n≔n+1;
9end


## 3.1. 1D Moving Mesh Partial Differential Equation

Letx 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 a,b and 0,1, respectively, is denoted by
(4)x=xξ,t,ξ∈0,1,with x0,t=a and x1,t=b. For a uniform mesh on the computational domain, given by ξj=j/N with j=0,1,⋯,N, a time-dependent mesh on the physical domain can be correspondingly obtained by setting xjt=xξj,t for all j. 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 xξ,t 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]:
(5)∂∂ξMxξ,t,t∂x∂ξ=0,with boundary conditions x0,t=a and x1,t=b. Here, Mx,t 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 ẋξ,t=∂/∂txξ,t. 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
(6)∂2ẋ∂ξ2=−1τ∂∂ξM∂x∂ξ,where τ is a positive parameter for adjusting the response time of mesh movement to the change of the monitor function Mx,t. With boundary conditions x0t=a and xNt=b, 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,
(7)xj+1n+1−2xjn+1+xj−1n+1−xj+1n−2xjn+xj−1nΔtn=−1τMj+1/2nxj+1n+1−xjn+1−Mj−1/2nxjn+1−xj−1n+1for j=1,2,⋯,N−1, where Δtn=tn+1−tn is the time step length, xjn≈xjtn is the numerical approximation of the jth mesh point at time tn, and Mj+1/2n=Mj+1n+Mjn/2 with Mjn=Mxjn,tn is the discrete monitor function on the jth mesh point at time tn. 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−Lx/2,Lx/2 and 0,1 be still denoted by x=xξ,t with x0,t=−Lx/2 and x1,t=Lx/2. Given a uniform mesh on the domain 0,1 with ξi=i/Nx for i=0,1,⋯,Nx, a time-dependent mesh on the domain −Lx/2,Lx/2 could be obtained by setting xit=xξi,t for all i. Similarly, by introducing a time-dependent one-to-one coordinate transformation y=yη,t with y0,t=−Ly/2 and y1,t=Ly/2 between 1D domains −Ly/2,Ly/2 and 0,1, a time-dependent mesh on the domain −Ly/2,Ly/2 could be obtained by yjt=yηj,t, where ηj=j/Ny with j=0,1,⋯,Ny is the uniform mesh on the domain 0,1. Then a time-dependent rectangular mesh on the physical domain Ω would be generated by setting the mesh point to be xit,yjt for all i and j.As stated in the previous subsection, bothxit and yjt can be determined from 1D MMPDE6 by utilizing appropriate monitor functions Mx,t and Gy,t, respectively, where Mx,t and Gy,t are functions of the 2D solution Tx,y,t, 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 transformationsx=xξ,t and y=yη,t between the physical coordinates x, y and the computational coordinates ξ, η, any function of x, y, and t can be expressed as a function in terms of ξ, η, and t, that is, fx,y,t=fxξ,t,yη,t,t. By the chain rule, it follows that
(8)∂f∂ξ=∂f∂x∂x∂ξ,∂f∂η=∂f∂y∂y∂η,∂f∂tξ,ηfixed=∂f∂tx,yfixed+∂f∂xẋξ,t+∂f∂yẏη,t.In order to distinguish the two partial derivatives with respect tot in the above expression, the notation ḟ, similar to the notation of mesh speed ẋ, is introduced for the first one, i.e., ḟ=∂f/∂tξ,ηfixed, and the other one is simplified to the original notation ∂f/∂t without causing confusion. Then, in the computational coordinates ξ, η∈0,1 and t>0, the original physical equation (1) becomes
(9)Ṫ−∂T∂ξ/∂x∂ξẋ−∂T∂η/∂y∂ηẏ−μ∂∂ξ∂T∂ξ/∂x∂ξ/∂x∂ξ+∂∂η∂T∂η/∂y∂η/∂y∂η=g~ξ,η,t,where μ=k/ρc is the thermal diffusivity and g~ξ,η,t=1/ρc∑l=1qglxξ,t,yη,t,t.The above equation can be discretized using the second-order central finite difference method on the uniform computational meshξi,ηj with i=0,1,⋯,Nx and j=0,1,⋯,Ny. This subsequently yields a system of ordinary differential equations of the form
(10)Ṫi,j−Ai,jtẋi−Bi,jtẏj−Li,jt=0,i=1,2,⋯,Nx−1,j=1,2,⋯,Ny−1,where Ti,j=Ti,jt=Tξi,ηj,t, xi=xit, yj=yjt, and
(11)Ai,jt=Ti+1,jt−Ti−1,jtxi+1t−xi−1t,Bi,jt=Ti,j+1t−Ti,j−1tyj+1t−yj−1t,Li,jt=μ2xi+1t−xi−1tTi+1,jt−Ti,jtxi+1t−xit−Ti,jt−Ti−1,jtxit−xi−1t+2yj+1t−yj−1tTi,j+1t−Ti,jtyj+1t−yjt−Ti,jt−Ti,j−1tyjt−yj−1t+g~ξi,ηj,t.Using the Crank-Nicolson method for temporal discretization, a full discretization, which has second-order time accuracy, can be obtained by(12)Ti,jn+1−Ti,jnΔtn−Ai,jn+Ai,jn+12xin+1−xinΔtn−Bi,jn+Bi,jn+12yjn+1−yjnΔtn−Li,jn+Li,jn+12=0,for i=1,2,⋯,Nx−1 and j=1,2,⋯,Ny−1. In the above equation, Ti,jn≈Ti,jtn is the numerical approximation of the temperature at ξi,ηj at time tn, equivalently at xin,yjn of the physical domain at time tn. As for Ai,jn, Bi,jn, and Li,jn, they are numerical approximations of Ai,jtn, Bi,jtn, and Li,jtn, respectively, and computed by substituting all time-dependent quantities with the corresponding numerical approximations in (11). Similarly, Ti,jn+1, Ai,jn+1, Bi,jn+1, and Li,jn+1 are corresponding numerical approximations at time tn+1.Supplemented with appropriate discretization of boundary condition (3), the linear system (12) can then be solved for all Ti,jn+1. Let us take the left boundary where x=−Lx/2 or equivalently ξ=0 as an example. If on the left boundary it is subjected to the Dirichlet boundary condition, we shall directly set
(13)T0,jn+1=T¯,for all j. Alternatively, if on the left boundary it is subjected to the Neumann boundary condition, which reduces to
(14)−k∂T∂x=q¯,we shall take the following discretization:
(15)T1,jn+1−T0,jn+1x1n+1−x0n+1+T2,jn+1−T0,jn+1x2n+1−x0n+1−T2,jn+1−T1,jn+1x2n+1−x1n+1=−q¯kfor all j, 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 xin+1 and yjn+1, 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 Min and Gjn.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 to0. In view of these, a linear combination of the arc-length monitor function and the curvature monitor function, which reads
(16)Mx,t=θ1+∂u∂x2+1−θ1+∂2u∂x224,is employed in our numerical experiments. Here, u=ux,t is a 1D function defined later by a certain average of the 2D temperature Tx,y,t with respect to y, and θ is the weight of the arc-length monitor function. Applying the central finite difference method to (16), one can obtain Min on xin by
(17)Min=θ1+ui+1n−ui−1nxi+1n−xi−1n2+1−θ1+2xi+1n−xi−1nui+1n−uinxi+1n−xin−uin−ui−1nxin−xi−1n24,where uin=uxin,tn. Apparently, it is enough to give uin in the computation of Min. Taking the whole 2D temperature field into consideration, the value of uin may be defined by
(18)uin=1Ny+1∑j=0NyTi,jn.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, Min is smoothed in our simulations by the strategy proposed in [28], i.e.,
(19)Min≔∑k=i−νi+νMkn2γ1+γ∣k−i∣/∑k=i−νi+νγ1+γ∣k−i∣,where γ>0 and ν≥0 are two smoothing parameters, given by γ=2 and ν=2 currently. Following the same approach, we can get Gjn by replacing i, xin, and uin in the right-hand side of (17), respectively, with j, yjn, and
(20)u~jn=1Nx+1∑i=0NxTi,jn.Then, Gjn would be smoothed with the same strategy of (19).Algorithm 1. Flowchart of the moving mesh algorithm for moving heat source problem.
Input: The end time tend, initial physical mesh (xi0, yj0) and initial temperature field Ti,j0.
Output: The final physical mesh (xin, yjn) and the corresponding temperature field Ti,jn.
1 Letn=0 and tn=0;
2while tn<tend do
3  Determine the time step∆tn;
4  Compute the 1D monitor functionsMin on xin for all i, and Gjn on yjn for all j, based on the current physical mesh (xin, yjn) and the corresponding temperature field Ti,jn;
5  Solve two linear systems of discretization of 1D MMPDE6 withMin and Gjn, respectively, to get two new 1D mesh xin+1 and yjn+1;
6  Construct the new physical mesh (xin+1, yjn+1);
7  Solve the system of discretization (12) to get the new temperature field Ti,jn+1;
8  Lettn+1=tn+∆tn and n≔n+1;
9end


## 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 Table1 are employed for the involved physical variables, and the plate is assumed to be homogeneous with the material density ρ=7.6×10−6, the heat capacity c=658, and thermal conductivity k=0.025. The room temperature 20 is adopted for both the initial temperature T0 and the boundary temperature T¯. When the Neumann boundary condition is taken into account, the boundary heat flux q¯ would be 0.001. Moreover, the time step length is given by Δtn=0.001, the parameter τ in MMPDE6 takes the value of 5×10−3, and the weight θ in the monitor function is set to be 0.05, if they are not explicitly pointed out below. For the rest of the parameters, they will be specified for each experiment individually.Table 1
Units of parameters for numerical experiments.
Length of the plate (Lx)
mm
Width of the plate (Ly)
mm
Time (t)
s
Temperature (T, T¯)
°C
Density (ρ)
Kg/mm3
Heat capacity (c)
J/kg/°C
Thermal conductivity (k)
W/mm/°C
Effective radius of heat source (rl)
mm
Maximum heat flux (Ql)
W/mm2
Boundary heat flux (q¯)
W/mm2


### 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 thex-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 Lx=100 and the width of Ly=50. 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 r1=2 and the maximum heat flux of Q1=5, is initially at the center of the right boundary and moves from right to left along x-axis with a constant speed of 2. It follows that α1t=50−2t and β1t=0.The test is simulated on the moving mesh with several different values ofNx and Ny. Numerical results are subsequently compared with the solutions obtained from the discretization (12) on the uniform mesh with various Nx and Ny. As presented in Figure 2, for the temperature profile along the heat source moving path, i.e., x-axis with y=0, at t=5, it can be seen that the solution on the moving mesh with Nx=50 and Ny=25 is much more accurate than the solution on the uniform mesh with the same Nx and Ny. In fact, it is comparable to the results in the uniform mesh with Nx=100 and Ny=50. In each time step, a single linear system of order Nx×Ny is required to be solved for the uniform mesh algorithm, whereas for the moving mesh algorithm, three linear systems of order Nx, Ny, and Nx×Ny, 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.Figure 2
Comparison of temperature profile (a) and their zooms (b, c) alongx-axis at t=5 obtained on the moving mesh with Nx=50, Ny=25 and the uniform mesh with various Nx and Ny.
(a)
(b)
(c)The transient 2D temperature field obtained by the moving mesh algorithm withNx=50 and Ny=25, together with the corresponding physical mesh, is presented in Figure 3, at time instances t=15 (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.Figure 3
The transient temperature field and the corresponding moving mesh att=15 (a), 25, 35, and 45 (d), respectively, for a single heat source moving along x-axis from right to left.
(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 Figure2. 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 lengthLx=Ly=100 is subjected to a single Gaussian point heat source, which moves along a circle of radius 15 with a constant speed in a counterclockwise direction. Specifically, the heat source has the effective radius of rl=2 and the maximum heat flux of Ql=15. Its moving path is set to be α1t=15cosπt/2 and β1t=15sinπt/2. Additionally, all boundaries of the plate are assumed to satisfy the Dirichlet boundary condition.This experiment is simulated by the moving mesh algorithm withNx=Ny=50 and the weight in the monitor function to be θ=0.2. The transient 2D temperature field at time instances t=1 (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.Figure 4
The transient temperature field and the corresponding moving mesh att=1 (a), 2, 3, and 4 (d), respectively, for a single heat source moving along a circle in a counterclockwise direction.
(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 ofNx and Ny. 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 Figure5, it would be stationary in the moving coordinate system that attaches to the moving heat source. Similar results can also be found in [29].Figure 5
The quasistationary temperature field for a single heat source moving along a circle in a counterclockwise direction.


### 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 alongx-axis in opposite directions, two heat sources moving along two intersecting straight lines, and three heat sources moving along three straight lines parallel to x-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 rl=2 and the maximum heat flux to be Ql=5, except for the last case where Ql=15.For the first case, the size of the plate is set to beLx=200 and Ly=100. The two heat sources are suddenly imposed on the position ±50,0, respectively, at the initial time, and then move along x-axis in opposite directions with a constant speed of 2. The resulting moving paths are α1t=−α2t=−50+2t and β1t=β2t=0.Obviously, the two heat sources will meet each other at time instancet=25. The simulation is performed by the proposed moving mesh algorithm with Nx=100 and Ny=50. The corresponding transient 2D temperature field as well as the physical mesh are presented in Figure 6, for time instances t=15 (a), 25, 35, and 45 (d), respectively. Additionally, the 1D temperature profiles along the heat source moving path at time instances t=15, 25, 35, 45, 55, and 65 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 1200. 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 650.Figure 6
The transient temperature field and the corresponding moving mesh att=15 (a), 25, 35, and 45 (d), respectively, for two heat sources moving along x-axis in opposite directions.
(a)
(b)
(c)
(d)Figure 7
The transient temperature profile alongx-axis at t=15 (a), 25 (b), 35 (c), 45 (d), 55 (e), and 65 (f), for two heat sources moving in opposite directions.
(a)
(b)
(c)
(d)
(e)
(f)For the second case, the plate is square with side lengthLx=Ly=100. The heat source moving paths are set to be α1t=−α2t=β1t=β2t=−25+2t. That is, the two heat sources are initially at the position ±25,−25, and move along the straight lines y=±x, respectively, with the constant speed of 2. Thus, they will meet each other at the original point 0,0 at time instance t≈17.678. The transient 2D temperature field and the corresponding physical mesh, obtained by the moving mesh algorithm with Nx=Ny=50, are plotted in Figure 8 for time instances t=10 (a), 17.678, 25, and 35 (d), respectively. Similar phenomena could be observed as the previous case.Figure 8
The transient temperature field and the corresponding moving mesh att=10 (a), 17.678, 25, and 35 (d), respectively, for two heat sources moving along the straight lines y=±x.
(a)
(b)
(c)
(d)For the last case, the plate is the same as the first case, i.e.,Lx=200 and Ly=100. The moving paths of the three heat sources are set to be α1t=α2t=α3t=100−20t, β1t=−β3t=20, and β2t=0, 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 20. The resulting transient 2D temperature field and the corresponding physical mesh by the moving mesh algorithm with Nx=100 and Ny=50 are shown in Figure 9 for time instances t=1 (a), 3, 5, 7, and 9 (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.Figure 9
The transient temperature field and the corresponding moving mesh att=1 (a), 3, 5, 7, and 9 (e), respectively, for three heat sources moving along the straight lines parallel to x-axis.
(a)
(b)
(c)
(d)
(e)


## 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 thex-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 Lx=100 and the width of Ly=50. 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 r1=2 and the maximum heat flux of Q1=5, is initially at the center of the right boundary and moves from right to left along x-axis with a constant speed of 2. It follows that α1t=50−2t and β1t=0.The test is simulated on the moving mesh with several different values ofNx and Ny. Numerical results are subsequently compared with the solutions obtained from the discretization (12) on the uniform mesh with various Nx and Ny. As presented in Figure 2, for the temperature profile along the heat source moving path, i.e., x-axis with y=0, at t=5, it can be seen that the solution on the moving mesh with Nx=50 and Ny=25 is much more accurate than the solution on the uniform mesh with the same Nx and Ny. In fact, it is comparable to the results in the uniform mesh with Nx=100 and Ny=50. In each time step, a single linear system of order Nx×Ny is required to be solved for the uniform mesh algorithm, whereas for the moving mesh algorithm, three linear systems of order Nx, Ny, and Nx×Ny, 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.Figure 2
Comparison of temperature profile (a) and their zooms (b, c) alongx-axis at t=5 obtained on the moving mesh with Nx=50, Ny=25 and the uniform mesh with various Nx and Ny.
(a)
(b)
(c)The transient 2D temperature field obtained by the moving mesh algorithm withNx=50 and Ny=25, together with the corresponding physical mesh, is presented in Figure 3, at time instances t=15 (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.Figure 3
The transient temperature field and the corresponding moving mesh att=15 (a), 25, 35, and 45 (d), respectively, for a single heat source moving along x-axis from right to left.
(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 Figure2. 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 lengthLx=Ly=100 is subjected to a single Gaussian point heat source, which moves along a circle of radius 15 with a constant speed in a counterclockwise direction. Specifically, the heat source has the effective radius of rl=2 and the maximum heat flux of Ql=15. Its moving path is set to be α1t=15cosπt/2 and β1t=15sinπt/2. Additionally, all boundaries of the plate are assumed to satisfy the Dirichlet boundary condition.This experiment is simulated by the moving mesh algorithm withNx=Ny=50 and the weight in the monitor function to be θ=0.2. The transient 2D temperature field at time instances t=1 (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.Figure 4
The transient temperature field and the corresponding moving mesh att=1 (a), 2, 3, and 4 (d), respectively, for a single heat source moving along a circle in a counterclockwise direction.
(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 ofNx and Ny. 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 Figure5, it would be stationary in the moving coordinate system that attaches to the moving heat source. Similar results can also be found in [29].Figure 5
The quasistationary temperature field for a single heat source moving along a circle in a counterclockwise direction.


## 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 alongx-axis in opposite directions, two heat sources moving along two intersecting straight lines, and three heat sources moving along three straight lines parallel to x-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 rl=2 and the maximum heat flux to be Ql=5, except for the last case where Ql=15.For the first case, the size of the plate is set to beLx=200 and Ly=100. The two heat sources are suddenly imposed on the position ±50,0, respectively, at the initial time, and then move along x-axis in opposite directions with a constant speed of 2. The resulting moving paths are α1t=−α2t=−50+2t and β1t=β2t=0.Obviously, the two heat sources will meet each other at time instancet=25. The simulation is performed by the proposed moving mesh algorithm with Nx=100 and Ny=50. The corresponding transient 2D temperature field as well as the physical mesh are presented in Figure 6, for time instances t=15 (a), 25, 35, and 45 (d), respectively. Additionally, the 1D temperature profiles along the heat source moving path at time instances t=15, 25, 35, 45, 55, and 65 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 1200. 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 650.Figure 6
The transient temperature field and the corresponding moving mesh att=15 (a), 25, 35, and 45 (d), respectively, for two heat sources moving along x-axis in opposite directions.
(a)
(b)
(c)
(d)Figure 7
The transient temperature profile alongx-axis at t=15 (a), 25 (b), 35 (c), 45 (d), 55 (e), and 65 (f), for two heat sources moving in opposite directions.
(a)
(b)
(c)
(d)
(e)
(f)For the second case, the plate is square with side lengthLx=Ly=100. The heat source moving paths are set to be α1t=−α2t=β1t=β2t=−25+2t. That is, the two heat sources are initially at the position ±25,−25, and move along the straight lines y=±x, respectively, with the constant speed of 2. Thus, they will meet each other at the original point 0,0 at time instance t≈17.678. The transient 2D temperature field and the corresponding physical mesh, obtained by the moving mesh algorithm with Nx=Ny=50, are plotted in Figure 8 for time instances t=10 (a), 17.678, 25, and 35 (d), respectively. Similar phenomena could be observed as the previous case.Figure 8
The transient temperature field and the corresponding moving mesh att=10 (a), 17.678, 25, and 35 (d), respectively, for two heat sources moving along the straight lines y=±x.
(a)
(b)
(c)
(d)For the last case, the plate is the same as the first case, i.e.,Lx=200 and Ly=100. The moving paths of the three heat sources are set to be α1t=α2t=α3t=100−20t, β1t=−β3t=20, and β2t=0, 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 20. The resulting transient 2D temperature field and the corresponding physical mesh by the moving mesh algorithm with Nx=100 and Ny=50 are shown in Figure 9 for time instances t=1 (a), 3, 5, 7, and 9 (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.Figure 9
The transient temperature field and the corresponding moving mesh att=1 (a), 3, 5, 7, and 9 (e), respectively, for three heat sources moving along the straight lines parallel to x-axis.
(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.



---

*Source: 6067854-2020-02-11.xml*