Numerical Investigation of Heat Conduction with Multiple Moving Heat Sources

This paper is concerned about the efficiently numerical simulation of heat conduction problems with multiple heat sources that are allowed to move with different speeds. Based on the dynamical domain decomposition upon the trajectories of moving sources, which are solved by a predictor–corrector algorithm, a non-overlapping domain-decomposed moving mesh method is developed. Such a method can not only generate the adaptive mesh efficiently by parallel computing, but also greatly simplify the discretization of the underlying equations without loss of accuracy. Numerical examples for various motions of sources are presented to illustrate the accuracy, the convergence rate and the efficiency of the proposed method. The dependence of the solution on the moving sources such as the types of motion and the distance between sources is numerically investigated. A blow-up phenomenon that occurs at multiple locations simultaneously can also be well observed for the case of symmetrically moving sources.


Introduction
The heat conduction transfer of problems with moving heat sources has attracted a great deal of interest by engineers in the past few decades.Typical applications include contact surfaces such as free-boundary solidification [1], and many metallurgical processes such as laser cutting and welding [2][3][4][5].Mathematically, this problem can be well modeled by the heat conduction equation with singular source terms, which utilize a time-dependent delta function to describe each highly localized and moving heat source (see [1,3,6,7] and references therein).It is well known nowadays that the solution of such a model equation is continuous and piecewise smooth with a jump for derivative of the solution when crossing each heat source [8].It has also been proved theoretically that blow-up may occur at some location in a finite time, if one of the source strengths can not be neglected and the corresponding source is traveling with a relatively low speed [7].Due to these singularities, numerical simulation to investigate the solution phenomenon, such as the prediction of the blow-up time, accurately and efficiently still remains big challenges.
When the singularity behaves in a time-dependent localized region, adaptive mesh methods are more popular than the numerical methods using a fixed mesh, which always leads to a huge computational cost to preserve accuracy, as the solution evolves singularity.Among the adaptive mesh methods, this paper is mainly concerned about the moving mesh method, which is easy to use and has been successfully applied to investigate the blow-up phenomenon [9,10] and other time-dependent localized singularities [11][12][13][14].The moving mesh, which continuously adjusts the locations of mesh points to resolve the solution as well as possible with a fixed total number of mesh points, can be generated adaptively by the moving mesh partial differential equations (MMPDEs) derived from the equidistribution principle [15,16].Based on the Schwarz domain decomposition method, one of the overlapping domain decomposition methods, the parallel generation of moving mesh has been studied in recent years [17][18][19].
For problems with moving heat sources, a second-order moving mesh method has been developed for one source case in [20], where the source trajectory is given explicitly and the delta function is approximated in the final discretization by appropriate information of jumps.Some modifications, to simplify the discretization without loss of accuracy, have been proposed and extended to multiple sources case in [21,22].However, all of these moving mesh methods were only presented for the one source case or the multiple sources case that all sources are moving with the same speed.Nevertheless, it becomes too complicated for [20,21] and even impossible for [22] to extend the corresponding method to the situation that the sources are moving with different speeds.Moreover, it is worth mentioning that these methods are usually failed to capture the blow-up phenomenon occurring at multiple locations simultaneously.Even worse, these methods can not be smoothly generalized to multi-dimensional situations, since it is difficult to define appropriate jumps of solution caused by the delta function as in a one-dimensional case.
Motivated by the above observation, this paper focuses on the development of the moving mesh method for problems with multiple moving heat sources, whose motions are controlled by different ordinary differential equations.The main idea and novelty is to decompose the domain dynamically according to the locations of moving heat sources, such that there exists a fixed mesh point moving exactly along the trajectory of each source, which is solved by a predictor-corrector algorithm.In each time-dependent subdomain, an MMPDE is applied to find the local equidistributed mesh.The model equation is then solved on the adaptive mesh assembling from local mesh on each subdomain.There are many advantages of the present method in comparison to the methods proposed in [20][21][22]:

•
Firstly, the adaptive mesh is optimal in the sense of equidistribution on each subdomain, and can be generated very efficiently, since the local equidistributed mesh on each subdomain can be obtained in parallel.

•
Secondly, it is found that the jump [ u] always vanishes on the specially designed adaptive mesh.Consequently, the discretization of the underlying model equation is further simplified while maintaining a second-order convergence rate in space.Moreover, the derived discretization is reinterpreted from a novel perspective, which shows some inspiration to extend the proposed method to a multi-dimensional case.

•
Thirdly, with the help of exact tracing of each heat source by a fixed mesh point, the blow-up phenomenon that occurs at multiple locations simultaneously has been successfully observed by the present numerical experiment of symmetrically moving sources.

•
Finally, various numerical experiments are carried out to investigate the solution behavior in terms of the motions of heat sources, such as the types of motion and the distance between the sources.The blow-up phenomena in terms of the distance of two moving sources are numerically discussed.It is found as expected that, if the sources are moving with high speeds, there may exist a critical distance, such that blow-up must occur if the distance of the sources is less than the critical distance, while blow-up is avoided when the distance is larger than this critical distance.
The rest of this paper is organized as follows.The model equation for multiple moving heat sources problem is briefly introduced in Section 2. The detailed numerical method is described in Section 3. Numerical examples are exhibited in Section 4 to validate the accuracy and the efficiency of the present method, and to investigate the solution behavior for various motions of heat sources.Some conclusions are given in the last section.

Model Equation
The temperature distribution in a reactive-diffusive material with multiple moving heat sources can be well described by the following equation [7,20]: with the initial and boundary conditions where u(x, t) is the temperature, with u t and u xx representing its first order partial derivative with respect to t and its second order partial derivative with respect to x, respectively.The function u 0 (x) is continuous with u 0 (x) → 0 as |x| → ∞, and q is the number of heat sources.The ith heat source, which is modeled by the Dirac delta function δ(x − α i (t)) with the source function F i (u), is moving along the curve α i (t), and its velocity is determined by where ψ i is a given function that may depend on the local temperature u(α i (t), t).Currently, the author focuses on the case that all sources do not intersect with each other during the time in consideration.Additionally, the source function F i (u) is usually a nonlinear function of u, and it has been proved [6,7] that blow-up may occur at some location in a finite time T, i.e., lim t→T max{u(x, t)} → +∞, if one of the source functions, namely F i 0 (u), satisfies and the corresponding source is moving with a sufficiently low speed, that is, |α i 0 (t)| ≤ κ, where κ is an appropriate constant.

Numerical Method
This section is going to develop a numerical method to solve the model Problems (1)-( 4) accurately and efficiently on a bounded domain [x l (t), x r (t)] that contains all heat sources.I will first introduce an adaptive mesh generated by the domain-decomposed moving mesh method and then give the detailed discretization of the underlying model equations on this special mesh.The final algorithm will then be presented.Without loss of generality, let us assume that x l (t)
In each physical subdomain [α i−1 (t), α i (t)], the moving mesh method is used to adaptively adjust the mesh such that the appropriately defined monitor function M = M(t, x, u) is equidistributed.The monitor function plays a key role in the moving mesh method [16], and will be given in Section 4 for our numerical experiments.Among the popular moving mesh partial differential equations (MMPDEs) proposed in [15] to equidistribute the monitor function, the so-called MMPDE6, which reads is employed as an example in this paper.Here, ẋ is the first order partial derivative of x with respect to t, and τ is a user-defined positive number for adjusting the time scale of the mesh equation.
With boundary conditions the MMPDE6 ( 5) can be solved by the following finite difference discretization is the average of discrete monitor function on the jth and (j + 1)th mesh points, and x n j is the numerical approximation of x(ξ j , t n ).When Equation ( 7) together with boundary Conditions (6) have been solved on all subdomains, we then get the new mesh G n+1 = {x n+1 j } for the whole domain [x l (t n+1 ), x r (t n+1 )] (see Figure 1 for the sketch of the above mesh generation).Again, it is not necessary to utilize the same monitor function and the same MMPDE in each subdomain.Thus, the present mesh generation strategy is very flexible.It is also pointed out that such an adaptive mesh can be obtained very efficiently by parallel computing based on domain decomposition methods [23].Moreover, it is optimal in the sense of equidistribution of the monitor function on each subdomain.Further advantages of such a mesh can be obtained in the discretization of the model Equation ( 1), which will be given in the next subsection.
Given the physical domain [x l , x r ], the physical discretization mesh G n and the corresponding solution u n j .
Compute the monitor function M.
Parallel Solve Equation ( 7) with boundary Conditions (6) on each subdomain Assemble the local equidistributed mesh of each subdomain to obtain the new mesh G n+1 .

Figure 1.
The sketch of domain-decomposed moving mesh generation.

Discretization on the Moving Mesh
Now, it is ready to introduce the discretization of the underlying Equations ( 1) and ( 4) on the adaptive moving mesh generated in the previous subsection.
For an arbitrary function f = f (x, t) = f (x(ξ, t), t), we have where f t and f x are the first order derivatives of f with respect to t and x, respectively.Consequently, the original Equation ( 1) can be rewritten into When x = α i (t), it is obvious that the above equation reduces to the following equation: Therefore, one may derive the discretization of Equation ( 8) based on the standard finite difference discretization of Equation ( 9), with modifications containing appropriate information of jumps when the mesh line of stencil mesh points intersects with one of the trajectories of heat sources.The jump of a given function ϕ along the smooth curve (θ(λ), ϑ(λ)) at the crossing point ( x, t) is defined by [20] [ϕ] ( x, t) := lim where (θ( λ), ϑ( λ)) = ( x, t).Then, from Equation ( 8), we have Physically, the value u(α i (t), t) changes smoothly as time evolves [20].It follows that the jump of the directional derivative of u(x, t) along the vector (α i (t), 1) is zero.By noting that, for the designed adaptive moving mesh, there is a mesh point that moves exactly along the trajectory of each heat source, we then obtain i = 0, 1, . . ., q − 1.Additionally, we can deduce that With the above information of jumps, the discretization of Equation ( 8), which has second-order convergence rate in space, is finally divided into two categories according to the position of the mesh point.Specifically, if the mesh point does not locate at the heat source, that is, j = 0, 1, . . ., N, but j = (i + 1)K, i = 0, 1, . . ., q − 1, we have the standard central difference discretization where x n j = x(ξ j , t n ), h n j = x n j − x n j−1 and u n j is the numerical approximation of the solution u(x n j , t n ).If the mesh point does exactly locate at the ith heat source, i = 0, • • • , q − 1, we have to incorporate the jump information into the discretization, which gives where j = (i + 1)K and ψ n+1 i is the numerical approximation of ψ i (t n+1 , α i (t n+1 ), u n+1 j ).For the detailed derivation of the above discretization, I refer the readers to [20].However, it is worth pointing out that, in [20], five cases must be considered for the discretization of Equation ( 8), even when there is only one heat source.
Remark 1.In the famous immerse boundary method [1,24], the delta function is first replaced by a suitable discrete approximation, which is usually the "hat function" with a fixed width support.Then, the resulting equation is discretized using standard numerical methods on the uniform grid.If now we approximate each delta function δ(x − α i (t)) in Equation (8) by the hat function defined as where h j (t) = x j (t) − x j−1 (t) and j = (i + 1)K, then applying the standard finite difference method for the resulting equation on the designed adaptive moving mesh G = {x j (t)} would also give the same discretization as Equations ( 13) and (14).
From the above point of view, first replacing the delta function by a mesh-dependent discrete approximation and then discretizing the resulting equation by standard numerical methods, one could obtain higher-order discretization both in space and time without much effort.It also shows some inspiration to extend the present method to multi-dimensional situation since the definition of jump is completely avoided, and this will be further studied in our future work.
Apparently, we need to know the source position α i (t n+1 ) before solving the system of Equations ( 13) and ( 14).If ψ i in Equation ( 4) is independent from u, one may only employ the following Crank-Nicolson scheme to update the source position α n+1 i ≈ α i (t n+1 ) as in [22,25].Otherwise, a predictor-corrector iteration will be employed to update the source position α n+1 i , the mesh points x n+1 j and the solution u n+1 j simultaneously.The details are given as follows.In the predictor step, let ψ n+1 i = ψ n i and solve Equation (15) to get an approximation α * i .Then, according to α * i , we solve Equations ( 7), ( 13) and ( 14) to obtain x * j and u * j .Finally, in the corrector step, we solve Equation (15) again with . This iteration is repeated until the convergence criterion is fulfilled.It is shown in Example 1 of Section 4 that such an algorithm is able to give satisfactory results even for the case that only one iteration is taken to update α n+1 i .To complete the problem, it remains to give the boundary conditions of the Equation (1) on the physical domain [x l (t), x r (t)].If the domain [x l (t), x r (t)] is sufficiently large, one would simply employ the zero Dirichlet boundary conditions, that is, u n+1 0 = u n+1 N = 0. Otherwise, the third order local absorbing boundary condition (LABC) proposed in [26], would be used.Here, s 0 is a user-defined parameter, the plus sign in "±" corresponds to the LABC at the right or left boundary.Using the coordinate transformation x(ξ, t), we can rewrite Equation ( 16) as Thus, we get the following central finite difference discretization: on the left boundary, and on the right boundary, where x −1 = 2x 0 − x 1 and x N+1 = 2x N − x N−1 .

Numerical Algorithm
This section is concluded by giving the whole algorithm in Figure 2 to solve the underlying Problems (1)-( 4).In the algorithm, we use t n+1 ≥ T as a stop criterion for a non-blow-up case and ∆t n < TOL for a blow-up case, where T is the termination time and TOL is a given small tolerance.
Given the initial mesh G 0 and initial solution u 0 .Let n = 0, t n = 0.
Determine ∆t n and let t n+1 = t n + ∆t n .
If ψ i depends on u?
Update the physical mesh G * (see Figure 1).

Numerical Examples and Discussion
Some numerical examples are presented in this section to verify the convergence rate, and to illustrate the accuracy and the efficiency of the proposed algorithm (see Figure 2).The solution behavior for various motions of heat sources are also investigated.
Example 1.We consider a nonlinear moving interface problem with the exact solution for some choice of ω 1 and ω 2 .The interface α 0 (t) is determined by solving the scalar equation so that u(x, t) is continuous across the interface.
Equation ( 21) has a unique solution [1] on the physical domain [0, 1] if ω 1 = 5π/4 and ω 2 = 7π/4, in which case, one has the initial interface position α 0 (0) ≈ 0.58333 and the motion of the interface satisfies the ordinary differential equation Based on the jump conditions, the source function F 0 (u(x, t)) is The monitor function used to generate the adaptive mesh in MMPDE6 ( 5) is taken as where θ = 0.5, ε = 10 3 /N 4 , and τ in MMPDE6 is set to be 10 −3 .In practical simulations, the monitor function is usually replaced by a smoother one [27] as where γ > 0 and ν ≥ 0 are two smoothing parameters, given by γ = 2 and ν = 2 in the following experiments.The uniform time step size is set to be ∆t n ≡ T L , and the termination time is taken as Our goal is to verify the present algorithm has a second-order convergence rate for space.Because the exact solution is known, we simply employ Dirichlet boundary conditions u(0, t) and u(1, t) at the boundaries.The numerical results with different numbers of N and L at the time T = 0.1 are listed in Table 1, where E d N is the maximum error of u obtained from [8], and the other errors are defined by in which u e and α e 0 are the exact solutions of Equations ( 20) and ( 21) respectively (α e 0 is obtained by the zero-finding MATLAB function fzero (MATLAB Release 2010a, The MathWorks, Inc., Natick, Massachusetts, MA, USA)).In addition, α0 is the solution of Equation ( 15) by one step of the predictor-corrector iteration.The numerical solutions U and Ũ are obtained by the algorithm shown in Figure 2 with α e 0 and α0 , respectively.The ratio is defined as the error in the coarse mesh divided by the one in the finer mesh.The first two columns in Table 1 show that α0 is a good approximation of the exact solution α e 0 , although we only adopt one step of the predictor-corrector iteration.Moreover, there is no remarkable difference between U and Ũ.Both of them are better than those obtained in [8].Finally, the numerical results illustrate that the convergence rate is second order in space.
Figure 3 shows the evolution of the mesh and the solution profiles for Problems (20) and ( 21) with N = 24.From the middle and the right subgraphs of Figure 3, it can be observed that the numerical solution (the solid line) are almost the same as the exact solution (the dots) at six different time stages t = 0, 0.02, 0.04, 0.06, 0.08, 0.1.The following examples are from moving source problems whose solutions may be blow-up [20][21][22].The physical domain is chosen as [−10, 10] with u(−10, t) = u(10, t) = 0 if it is not specifically pointed out.The initial value is given by and the source functions F i (u(x, t)) are specified by The monitor function in MMPDE6 (5) reads where ε ∈ (0, 1], p > 0, θ i ∈ (0, 1) and , where µ = 10 −3 .Additionally, the Newton iteration method is adopted to solve the resulting nonlinear system with the tolerance tol = 10 −8 .
Example 2 (Linear moving sources).We consider that all sources move with a constant velocity k, Since the proposed algorithm is trivial for the problem with multiple sources, only numerical results for q = 1, 2 are presented here.In these cases, the solution phenomenon with k = 1, 2 has been investigated in [20][21][22].Our numerical results give a satisfactory agreement with those shown in [20][21][22].In more detail, the blow-up occurs with k = 1, but it is avoided with k = 2 in the case of q = 1.The evolution of the moving mesh and the solution profiles in the case of k = 2 and q = 1 can be found in Figure 4, where the simulation parameters are given by τ = 10 −3 , θ 0 = 0.9, θ 1 = 0.1, and ε = 10 3 /N 4 .Now let us add another heat source with d 1 = 2.5 and k = 2; then, the phenomenon of blow-up can be observed again.With the simulation parameters τ = 5 × 10 −4 , θ 0 = θ 1 = 0.3, θ 3 = 0.4, p = 2, and ε = 10 −5 , it is shown from Figure 5 that the blow-up occurs at the position of the first source, i.e., x = 4.079417297361286, with the maximum value u max = 3.16 × 10 6 , and the corresponding blow-up time is 2.039708648680643.This result is very close to the one obtained in [22], where the blow-up time for this case is 2.039311846765550.
The blow-up phenomena of the solutions in terms of d 1 for the two sources case are also investigated.Since the blow-up occurs in two sources case while it is avoided in one source case when k = 2, it is expected that the blow-up time will be delayed as the distance d 1 increases, and there exists a critical distance, such that the blow-up does certainly occur if d 1 is less than the critical distance, while the blow-up will be avoided if d 1 is larger than the critical distance.As shown in Table 2, the blow-up time is delayed as expected as d 1 increases, yet the blow-up will still be observed even for d 1 = 100 when k = 2.However, it can be seen that, in the case d 1 ≥ 0.85 for k = 3 and d 1 ≥ 0.04 for k = 4, the blow-up is avoided and the maximum of the solution u max would converge to the value listed in Table 2 as time evolves.A further computation using N = 1500 shows that the critical distance for k = 4 is between 0.03902 and 0.03903.Example 3 (Sin-type moving sources).We now consider the case of two sources, in which the sources move periodically with the same speed For periodically moving sources, only one source case has been studied in [21].Our numerical results are shown in Figure 6 for A = π.The simulation parameters are the same as those in Example 2. It is shown that the blow-up occurs at the location of the second source with the maximum value u max = 3.16 × 10 6 , and the corresponding blow-up time is 1.689611393639939.This result is also very close to the one obtained in [22], where the blow-up time for this case is 1.688671486061044.Example 4 (Symmetrically periodically moving sources).We consider the case of two sources, which move periodically and symmetrically.The motions are described by where A = π.
To the best of my knowledge, there are no theoretical and numerical results for multiple sources with different speeds at present.For the two sources moving periodically and symmetrically, it is anticipated that if the blow-up occurs, then it occurs at the locations of both sources at the same time.However, it is difficult to simulate this symmetrical blow-up phenomenon due to numerical perturbation.Indeed, the globally moving mesh strategies proposed in [20,21] usually fail to capture the symmetrical blow-up phenomenon.However, with the benefit of the domain-decomposed adaptive mesh generation, one can simulate this symmetrical blow-up phenomenon very well.As shown in Figure 7, the blow-up occurs as expected at the locations of both sources in the time 2.496881990359248, where the maximum value of the solution is u max = 3.16 × 10 6 .For a comparison, we also present the simulation results with LABCs (17) on a much smaller physical domain [α 0 (t) − 4.0, α 1 (t) + 4.0] in Figure 8.It can be seen that the blow-up time 2.496370241342059 is very close to that in Figure 7, which shows the effectiveness of the LABCs.That is, a smaller physical domain with LABCs, indicating smaller computational cost, can be employed in the numerical investigation of the moving heat sources' problems.

Conclusions
A domain-decomposed moving mesh method has been developed for the heat conduction problem with multiple localized sources, which can be moving with different speeds.The special adaptive moving mesh can be generated very efficiently by parallel computing, and the discretization of the model equation on it is simplified a lot without loss of accuracy.A predictor-corrector algorithm is proposed to capture the trajectories of the heat sources and the solution of the model equation accurately.To the best of my knowledge, this is the first time to investigate the solution phenomenon for the problem of multiple heat sources with different speeds.With the proposed method, we have investigated the dependence of the solution on the moving sources such as the types of motion and the distance between sources.The symmetrical blow-up phenomenon in the case that two sources are moving symmetrically and periodically is also successfully simulated.
Finally, the derived discretization of the underlying equation is reinterpreted from a new point of view, which shows that the same discretization can also be derived by first replacing the delta function with a mesh-dependent hat function, and then applying standard numerical methods to the resulting equation.Since the jumps of solution are completely avoided in the new approach, one may extend the present method to multi-dimensional problems with moving localized heat sources.Actually, this work is ongoing and would be reported elsewhere soon.For the case that the moving localized heat sources may be crossing with other sources, it is obvious that the domain-decomposed moving mesh method could not be applied directly.However, the idea, approximating the delta function by a mesh-dependent hat function, may still work and give an accurate discretization, which will be considered in my future work.

Figure 3 .
Figure 3. Mesh trajectories (left), the profiles of u in terms of x (middle) and ξ (right) with N = 24.The solid lines are the numerical solution and the dots are the exact solution.

Figure 4 .
Figure 4. Mesh trajectories and the profiles of u from t = 0 to t = 1.0 for one source case with N = 100.

Figure 5 .
Figure 5. Mesh trajectories and the profiles of u with N = 150 in the case of two sources.

Figure 6 .
Figure 6.Mesh trajectories and the profiles of u for q = 2, A = π with N = 150.

Figure 7 .
Figure 7. Numerical results for symmetrically periodically moving sources with N = 150.