Semi-analytic approximation of the temperature field resulting from moving heat loads

Moving heat load problems appear in many manufacturing processes, such as lithography, welding, grinding, and additive manufacturing. The simulation of moving heat load problems by the finiteelement method poses several numerical challenges, which may lead to time consuming computations. In this paper, we propose a 2D semi-analytic model in which the problem in two spatial dimensions is decoupled into three problems in one spatial dimension. This decoupling significantly reduces the computational time, but also introduces an additional error. The method is applied to a wafer heating example, in which the computational time is reduced by a factor 10 at the cost of a 4% error in the temperature field. 2018 The Authors. Published by Elsevier Ltd. This is an openaccess article under the CCBY license (http:// creativecommons.org/licenses/by/4.0/).

The basis of the theory for moving heat sources was developed by Rosenthal [1,22], who observed that when the path of the heat load is long enough, the temperature distribution around the source soon becomes constant. Assuming constant material properties, Rosenthal developed closed-form analytic expressions for these quasi-stationary temperature fields resulting from point, line, and plane heat sources. Although Rosenthal's analysis provides valuable estimates, transient effects and position or temperature-dependent coefficients are important in many applications. In these situations the problem is solved by Finite Element (FE) analysis (see for example [2,3,16]).
Solving a moving heat load problem by the FE method poses several numerical challenges. One problem is that by fixing the coordinate frame to the heat load we obtain a convectiondiffusion problem. It is well known that the FE discretization of such problems may result in spurious oscilations [23]. Spurious oscilations can be prevented in two ways. In the first approach, the mesh size in the direction of the velocity of the moving load is chosen smaller than 2D=v, where D [m 2 /s] denotes the thermal diffusivity of the material and v [m/s] denotes the velocity of the moving load [23]. Note that this approach is computationally demanding when the velocity v is high. In the second approach, upwinding schemes [23,24] are used. These schemes prevent spurious oscilations at the cost of an increased discretization error.
Another problem is that the area in which the heat load is applied is typically small. This makes both the spatial and temporal discretization of such problems computationally demanding. For example, Zhang et al. [15] report that for a Gaussian heat distribution, the mesh size should be at least twice as small as the radius of the heat distribution and at least two time steps are needed for the time that the heat load travels along one element.
Because of these considerations, many problems require a small mesh size. For a static mesh, this mesh size needs to be used in the whole region through which the heat load travels, which results in models with many Degrees of Freedom (DOFs). To  proposed [3,4,17], which lead to significant reduction in computational effort. Note that these schemes require some cost for updating the mesh and that the temporal discretization remains challenging, since the adaptive mesh will keep the mesh size near the source small.
For problems with constant coefficients, spatial discretization can be avoided by semi-analytic methods [5][6][7][8]25]. In these methods, the temperature field is expressed as the convolution of the fundamental solution of the heat equation and the applied heat load. The convolution over space can typically be solved analytically, so that only numerical evaluation of the convolution over time remains. This is still a computationally intensive operation when the solution is evaluated on a fine grid.
In this paper, we propose a novel semi-analytic approximation method to reduce the computational cost of 2D transient moving load problems with constant coefficients. We construct a semianalytic approximation in which we decouple the problem in two spatial dimensions into three problems in one spatial dimension. This significantly reduces the computational cost, especially on fine grids. The proposed method is demonstrated by an example from precision engineering, more specifically for a wafer heating problem.
The remainder of this paper is structured as follows. In Section 2, the semi-analytic approximation for the temperature field is introduced on an infinite domain. In Section 3, we give a physical interpretation of the semi-analytic approximation. In Section 4, we discuss the modeling of edge effects and repetitive scanning patterns, which are typically encountered in lithography and additive manufacturing. In Section 5, we apply the developed techniques in a wafer heating example. In Section 6, the conclusions are presented and the results are discussed.

Problem formulation
We consider heat conduction in a thin infinite plate with thickness H [m] and constant material properties (see Fig. 1 [W/m 2 K], respectively. Because the plate is thin, the temperature gradient along the thickness of the plate can be neglected. The resulting temperature field T 2D ¼ T 2D ðx; y; tÞ [K] relative to a reference temperature T r satisfies the heat equation, see for example [22] qcH  Fig. 1 can be written in this form by taking block functions for XðxÞ and YðyÞ. Such a heat load has been considered in laser hardening [14] and will also be considered in the lithography example in Section 5. Also the Gaussian heat distribution considered in many applications (see for example [2,3,15,18]) is of the form in (2). Observe that Q moves with a constant velocity v in positive ydirection.

The approximate solution
We will introduce an approximate solution by simplifying the integral in (8). Note that the only factor in (8) that depends on x is Nðx; sÞ. In our semi-analytic approximation, we move this factor outside the integral. To this end, we approximate Nðx; sÞ by a first-order Taylor series expansion around s ¼ t Ã , where t Ã does not depend on s and is not determined yet Nðx; sÞ ¼ Nðx; t Ã þ ðs À t Ã ÞÞ % Nðx; t Ã Þ þ ðs À t Ã Þ @N @s ðx;sÞ¼ðx;t Ã Þ : ð11Þ When we substitute this approximation back into (8), we obtain a semi-analytic approximationT 2D of the true temperature field T 2D as The error in the semi-analytic approximationT 2D originates from the Taylor series approximation of Nðx; sÞ in (11), which can be expressed as (see e.g. [27]) wheret Ã is an (unknown) point between t Ã and s. Assuming the difference between t Ã and s is small, we may approximate the error iñ t Ã by t Ã , so that (15) becomes the first higher-order term that is omitted in the Taylor expansion (11). Based on (15) witht Ã replaced by t Ã , the error in the semianalytic approximation is estimated as We choose t Ã to minimize R 1D , so we compute where T 1D and T 1D are as in (13) and (14), respectively. We thus see that we should choose t Ã as Note that this choice of t Ã cancels the second term on the Right Hand Side (RHS) of (12), so that the semi-analytic approximation becomes T 2D ðx; y; tÞ ¼ T 1D ðy; tÞNðx; t Ã ðy; tÞÞ: We note that the semi-analytic approximationT 2D is welldefined with t Ã defined by (19). To see this, recall that Y and H are nonnegative, so that (9) shows that f is nonnegative. Hence 1D in (14) shows that When T 1D ðy; tÞ > 0, we may divide the above inequalities by T 1D ðy; tÞ, from which we find 0 6 t Ã ðy; tÞ 6 t: ð22Þ In case T 1D ðy; tÞ ¼ 0, (21) shows that T ð1Þ 1D ðy; tÞ ¼ 0. Now (12) shows thatT 2D ðx; y; tÞ ¼ 0 for any choice of t Ã . So when the quotient in (19) is not well-defined, the value of t Ã does not influence the resulting semi-analytic approximationT 2D .
Since the integral expressions in (10), (13), and (14) can generally not be solved explicitly, it will be convenient to express Nðx; sÞ; T 1D ðy; tÞ, and T ð1Þ 1D ðy; tÞ as the solutions of Partial Differential Equations (PDEs) in one spatial dimension. Using that Uðx; sÞ is the fundamental solution for the heat equation in one spatial dimension, we see that Nðx; sÞ in (10) is in fact the solution to the initial value problem and by looking back at (9) that T 1D in (13) is the solution to the PDE with zero initial conditions. To see how we can compute T ð1Þ 1D ðy; tÞ as the solution of a PDE, note that (14) can be rewritten as The integral in (26) can be interpreted as the solution to the PDE @T ð1cÞ 1D  (20).
Note that this procedure requires to solve three PDEs in one spatial dimension ((23), (24) and (25)), whereas (3) is a PDE in two spatial dimensions. This means the procedure to computẽ T 2D will be much more efficient than computing T 2D by discretizing (3). For example, when we consider the rectangular grid in Fig. 2 with N x grid points in the x-direction and N y gridpoints in the ydirection, the spatial discretization of (3) yields N x N y Ordinary Differential Equations (ODEs) whereas the spatial discretization of (24) and (23) only leads to N x þ 2N y ODEs.
In some cases, a closed form analytic expression for the solution for Nðx; sÞ can be obtained by computing the integral in (10) directly. A closed form analytic expression for T 1D ðy; tÞ and T ð1Þ 1D ðy; tÞ in (13) and (14) is impossible to obtain for most practical situations.
It is important to note that, at a certain level of accuracy, the truncation of the Taylor series in (11) introduces an error that cannot be reduced further by refining the 1D mesh. To reduce this error, more terms in the Taylor series expansion (11) should be considered, which leads to higher-order approximations. Because the semi-analytic approximationT 2D in (20) is based on the firstorder Taylor series approximation in (11),T 2D in (20) will be called the first-order semi-analytic approximation. The nth-order semianalytic approximation follows by substitution of an nth-order Taylor series approximation for Nðx; sÞ around s ¼ t Ã in (8). In particular, the second-order semi-analytic approximation is found by adding (16) to the first-order semi-analytic approximationT 2D in (20). The computation of an n-th order (n P 2) semi-analytic approximation is addressed in Appendix A.

Physical properties and interpretation
The (first-order) semi-analytic approximationT 2D defined in the previous section is build from three functions: T 1D ; N, and t Ã . These three functions have a clear physical interpretation in terms of the original problem.
For the physical interpretation, note that integrating Nðx; sÞ in where we applied a change of variables x 00 ¼ x À x 0 and changed the order of integration. Since the integral of Uðx; sÞ from x ¼ À1 to for all s P 0. Integrating the formula for T 2D in (8) over x and using (13) and (29), we see that By integrating (20) over x using (29), we also find Since the RHS of (30) is equal to the RHS of (31), we see that T 1D is chosen such that the internal energy of the semi-analytic approximation on lines in x-direction qcH R þ1 À1T 2D dx is matched with the internal energy of the exact solution along lines in the x-direction qcH R þ1 À1 T 2D dx. In particular, (29) shows that the shape functions Nðx; sÞ describing the dependence on x are normalized such that the integral over x is constant.
We can find a physical interpretation of N and t Ã when we assume that Yðy À vtÞ ¼ dðy À vtÞ is a Dirac delta and HðtÞ 1. In this case the heat load is applied on a line in the x-direction (see Fig. 3) and does not vary over time.
To find a physical interpretation for N, we think about which profiles we would expect to see in the x-direction on a line y ¼ y 0 . As observed by Goldak et al. [2], the heat transport into the y-direction may be neglected when the velocity v is high. Looking back at the original problem (3), this means that we neglect the term D@ 2 T 2D =@y 2 . On a line y ¼ y 0 , we thus expect that the profiles N y 0 ðx; tÞ in the x-direction are solutions of For zero initial conditions at t ! À1, the solution to this equation is N y 0 ðx; tÞ ¼ e hðy 0 =vÀtÞ 0 for t < y 0 =v where N is the solution of (23). The factor e hðy 0 =vÀtÞ does not depend on x, so that it does not change the shape profiles that are observed in the x-direction and therefore the profiles in the x-direction we expect to see are indeed snapshots from Nðx; sÞ.
For the interpretation of t Ã ¼ t Ã ðy; tÞ, it will be convenient to introduce a coordinate frame fixed to the heat load defined by ðx; f; tÞ ¼ ðx; y À vt; tÞ: Now observe that the semi-analytic approximation in (20) uses Nðx; t Ã Þ as the profile in the x-direction. Furthermore, the expected profile in (33) is Nðx; t À y 0 =vÞ. We thus expect t Ã to be where f 0 ¼ y 0 À vt. The expression for t Ã phys in (35) in fact describes the time since the heat load has arrived at y ¼ y 0 . This can be seen from Fig. 3: since the heat load is located at y ¼ vt (because Yðy À vtÞ ¼ dðy À vtÞ), the time the heat load arrives at y 0 is t 0 ¼ y 0 =v. At time t, the time since the heat load has arrived at y ¼ y 0 is t À t 0 ¼ t À y 0 =v, which is precisely t Ã phys . This interpretation is indeed closely related to the formula for t Ã in (19). For HðtÞ 1 and Yðy À vtÞ ¼ dðy À vtÞ, it easy to see from (9) that f ðy; t; sÞ ¼ f ðf þ vt; t; sÞ is only a function of f and s. Using this observation we can then compute lim t!1 t Ã ðf þ vt; tÞ ¼ where the integrals have been computed directly in MAPLE [28]. By expanding (36) in a first order Taylor expansion in 1=v around 1=v ¼ 0, we see that the expression in (36) approaches jfj=jvj for v ! 1. Since the choice of the line y ¼ y 0 was arbitrary in the derivation of t Ã phys , we may write y instead of y 0 in (35) and we find that t Ã phys ¼ Àf=v ¼ jfj=v for f < 0. We conclude that t Ã may be interpreted as the time since the heat load has arrived when Yðy À vtÞ ¼ dðy À vtÞ; HðtÞ 1, and v is large.

Edge effects
The results so far have been derived for a heat source moving in an infinite domain. However, in a practical situation such as the example in Section 5 we typically have a bounded domain.
Assuming the solution T R 2 ¼ T R 2 ðx; y; tÞ (for which we can compute a semi-analytic approximation using the approach from Section 2) on the infinite domain ðx; yÞ 2 R 2 is available, we discuss how we can construct solutions on the bounded domain.
The effect of one straight edge with a perfectly insulated or zero temperature boundary condition can be included by the method of images [1,29]. Fig. 4 (@T=@x ¼ 0, i.e. perfectly insulated edge) and For the wafer heating application, we are interested in using these reflection type of arguments for a circular domain with radius R and with perfectly insulated edges. To this end T R 2 is written in cylindrical ðr; hÞ-coordinates defined by x ¼ r cosðhÞ and y ¼ r sinðhÞ. To approximate the solution on the bounded domain r 6 R, we propose to reflect in the edge r ¼ R along radial lines according to the formula It is important to note that (37) does not give the exact solution on the bounded domain r 6 R, in contrast to the reflection along a straight edge. To make this more precise, note that for T in (37) to give the exact solution we need that (1) T should satisfy the boundary condition and (2) the reflected part T R 2 ðR 2 =r; h; tÞ should be a solution of the heat Eq. (3) with zero source term. To check the boundary condition, we compute @T @r so that T in (37) indeed satisfies the perfectly insulated boundary condition. However, the reflected partTðr; h; tÞ ¼ T R 2 ðR 2 =r; h; tÞ is not a solution of the homogeneous heat equation Note, however, that close to the edge r ¼ R the factor r 4 =R 4 % 1, meaning thatT is 'almost' a solution of the heat equation (39) near the edge. Since T R 2 ! 0 for r ! 1, we see that the reflected part T ! 0 when r ! 0, meaning thatT will trivially satisfy (39) in the center of the disc r ¼ 0. Combining these two observations we see that the reflection formula is accurate whenT is close to zero in the region where r 4 =R 4 is significantly different from 1.

Repetitive scanning patterns
The fact that we have been considering solutions on the infinite domain ðx; yÞ 2 R 2 is particularly useful when we are dealing with repetitive patterns, which occur in many applications such as additive manufacturing and lithography. As an example consider the expose pattern in Fig. 6, where the same heat load is applied to each of the 64 rectangles (fields), repeated in different locations in space and time. This means that the applied heat load can be written in the form Qðx; y; tÞ ¼ X Nr i¼1 Q 0 ðw i ðx; y À vðt À s i ÞÞ; t À s i Þ; where Q 0 ðx; y; tÞ denotes the heat load applied in one field which is repeated N r times, w i : R 2 ! R 2 describes the translation and rotation in space, and s i describes a shift in time. Since (2) on the infinite domain is linear and rotation, translation, and time invariant,   the temperature response resulting from the heat load in (41) can be obtained from the superposition principle as Tðx; y; tÞ ¼ X Nr i¼1 T 0 ðw i ðx; y À vðt À t i ÞÞ; t À t i Þ; where T 0 ðx; y; tÞ is the response resulting from Q 0 ðx; y; tÞ. This formula shows that we can cheaply obtain the solution T on the infinite domain, since we only need to compute the response T 0 for one field. We then obtain the solution on the bounded domain by applying the reflection arguments from Section 4.1.
We thus propose the following computational procedure to compute the temperature response for repetitive scanning patterns: This procedure is computationally much more efficient than discretizing (1) for the complete heat load Q. The most important reason is that the area through which the heat load Q 0 travels is much smaller than the area through which the complete heat load Q travels. This means that the computation of T 0 requires a fine mesh in a much smaller area than would be required to compute the response for the complete heat load Q. Another advantage is that Q 0 is applied during a shorter period of time than Q, so that the cost for the temporal discretization is also reduced. Finally, if we use the semi-analytic approximation discussed in Section 2, for the heat load applied in one field only N x þ 2N y ODEs need to be discretized instead of N x N y ODEs.

Example
We illustrate our method for a lithography example from the semiconductor industry. In this example, we consider the examplatory expose pattern in Fig. 6. Each of the 64 fields in the figure (fields) are scanned consecutively by a uniform heat load applied in a smaller rectangular area (the slit). The scanning of one field is shown in Fig. 7. We use the method from Section 4 to obtain the complete temperature field, so we will first focus on the computation of the response for one field T 0 , which we will then use to assemble the complete temperature field T. The considered parameter values used to generate the results are given in Table 1.
To verify the accuracy of the semi-analytic approximation, we compare the solution T 0 to FE solutions on several grids. We discretize (1) using rectangular 4-node bilinear quadrilateral elements. Around the area in which the heat load is applied we use square elements with length L e [m] and farther away from the heat load we increase the element size to efficiently approximate the solution on the infinite domain. Starting from the coarsest mesh considered for L e ¼ 6:4 mm shown in Fig. 7, we subdivide each element into four smaller ones until we reach L e ¼ 0:1 mm. For an FE simulation with L e ¼ 0:1 mm on the domain in Fig. 7 with perfectly insulated edges, the temperature increase on the edges is below 0.33% of the maximal temperature that occurs during the simulation, indicating that solutions computed on this domain will closely resemble the infinite domain solutions. The heat load is applied between t ¼ 0 and t ¼ t 1 ¼ 0:115 s and the time in the simulation runs from t ¼ 0 to t ¼ 3 s, so that also the passive cooling that occurs after the heat load has been applied is included in the simulation. The code for the spatial discretization is written in MATLAB, and the time integration is done using MATLAB's ODE solver ode15s with the default tolerances (a relative tolerance of 10 À3 and absolute tolerance of 10 À6 ).
The semi-analytic approximationT 2D is computed on the same grids as the FE solutions. The heat load in Fig. 7 can be written in the form (4): we set XðxÞ equal to 1=L for x 2 ½ÀL=2; L=2 and zero otherwise, YðfÞ equal to 1=W for f 2 ½ÀW=2; W=2 and zero otherwise, and Q ðtÞ is equal to 1 [W] for 0 6 t 6 t 1 and zero afterwards. For the considered XðxÞ, we can determine a closed-form analytic expression for Nðx; t Ã Þ from (10), resulting in where erf denotes the error function. To compute T 1D ðy; tÞ and T 1D ðy; tÞ we discretize (24) and (27), where we use the same grid in the y-direction as for the 2D FE solution. The spatial discretization is programmed in MATLAB and uses 2-node linear elements and the temporal discretization is again done by MATLAB's ode15s with the default tolerances. It is important to note that for larger grid sizes the discretization introduces spurious oscilations (see [23]) in T 1D ðy; tÞ and T ð1cÞ 1D ðy; tÞ, which can result in negative values of t Ã . Since (43) is only defined for t Ã > 0, we set Nðx; t Ã Þ ¼ Nðx; 0Þ ¼ XðxÞ when negative values of t Ã appear. Since Fig. 7. The scanning of one field on the wafer with the coarsest mesh that is considered. negative values of t Ã only appear when T 1D is small, (20) shows that this does not affect the semi-analytic approximation much. Fig. 8 shows the relative error in the temperature field for the FE solution and the first-and second-order semi-analytic approximation. The FE solution with L e ¼ 0:1 mm is used as reference. The relative error is computed by taking the maximum over space and time of the absolute value of the error, divided by the maximum of the reference solution over space and time. The error in FE solution decreases at a constant rate, but the error in the first-order semi-analytic approximation stops decreasing at L e ¼ 1:6 mm. At this point, the error in the first-order semi-analytic approximation is no longer dominated by discretization errors but by the error introduced by the truncation of the Taylor series approximation in (11). This limits the accuracy of the first-order semi-analytic approximation to about 4%. For the second-order semi-analytic approximation, the error stops decreasing at L e ¼ 0:4 mm and the accuracy is limited to about 1.5%. Fig. 9 shows the snapshot of the error profile of the first-order semi-analytic approximation at the moment t ¼ t m ¼ 0:1036 s at which the maximal error is observed. The relative error in Fig. 9 is computed by dividing the observed error by the maximal temperature observed during the simulation for the FE model with element size L e ¼ 0:1 mm. We see that the maximal relative error is indeed below 4% and that the error occurs near the left and right side of the field, so near x ¼ AEL=2. The location of the maximal error is in agreement with the leading term of the error in the first-order semi-analytic approximation in (16), which contains the factor @ 2 N @s 2 ðx; t Ã Þ. Since Nðx; sÞ is the solution of (23), we indeed expect that @ 2 N @s 2 ðx; t Ã Þ is the largest near discontinuities in the initial condition Nðx; 0Þ ¼ XðxÞ, which are indeed located at x ¼ AEL=2. Note that for smoother initial conditions XðxÞ, such as a Gaussian heat distribution, smaller errors as well as faster convergence may be expected when adding higher order terms.
We expect that the second-order semi-analytic approximation, which is obtained by adding (16) to the first-order semi-analytic approximation, will be more accurate, which is confirmed by Fig. 8. Since the coefficient R 1D ðy; t; t Ã Þ is nonnegative because of (17), the expected error profile in the x-direction for the firstorder semi-analytic approximation will be proportional to À @ 2 N @s 2 ðx; t Ã ðy; tÞÞ. Fig. 10a now shows the observed error profile in the first-order semi-analytic approximation in the x-directioñ T 2D ðx; y m ; t m Þ À T 2D ðx; y m ; t m Þ and Fig. 10b shows the expected error profile in the x-direction À @ 2 N @s 2 ðx; t Ã ðy m ; t m ÞÞ, where y m ¼ 32 mm and t m ¼ 0:1036 are such that maximal error occurs at the line y ¼ y m at time t ¼ t m . Indeed, the shape of both profiles is very similar. This suggests that the expression in (16) gives a good indication of the observed error. Fig. 8 confirms that when we compute the second-order approximation by adding the term R 1D @ 2 N @s 2 in (16) to the first-order approximation, the relative error is reduced to 1.5%. Fig. 11 shows the CPU running times needed to obtain the temperature field for the FE solution and the first-and second-order semi-analytic approximations. At L e ¼ 1:6 mm, the first-order semi-analytic approximation is computed almost 10 times faster than the FE solution, while a similar accuracy is obtained, see Fig. 8. Computing the second-order semi-analytic approximation increases the computational cost because one additional PDE in one spatial dimension must be solved and because the time integration using MATLAB's ode15s now requires stricter tolerances    10. The error in the first-order semi-analytic approximation (computed for an element size of Le ¼ 0:1 mm) along the line y ¼ y m ¼ 32 mm at the moment t ¼ tm ¼ 0:1036 s compared to À @ 2 N @s 2 ðx; t Ã ðy m ; tmÞÞ, which is the situation where the maximal error in the semi-analytic approximation is observed.
(the absolute and relative tolerance were set to 10 À8 for the second-order semi-analytic approximation). Nevertheless, at L e ¼ 0:8 mm the second-order semi-analytic approximation is still computed 10 times faster than the FE solution while a comparable accuracy is obtained, see Fig. 8.
From the response of one field we can now assemble the response for the complete scanning pattern using the approach from Section 4. A snapshot of the resulting temperature field is shown in Fig. 12 at the moment that the scanning of field nr. 35 (see Fig. 6 for the field numbering) is almost complete. The temperature in Fig. 12 is normalized w.r.t. the maximal temperature observed at t ¼ 4 s. The computational advantage of the method from Section 4.2 is clear in this case: computing the temperature field on a 0:5 Â 0:5 mm 2 mesh is practically impossible in a direct FE model, because the computation time exceeds several days, whereas the method from Section 4.2 only takes less than 200 s to compute the solution in 2000 time steps.

Conclusions and discussions
We have introduced a semi-analytic approximation for the calculation of the 2D temperature field resulting from a moving heat load, by decoupling the problem in two spatial dimensions into three problems in one spatial dimension. Especially on fine meshes, this leads to a significant reduction in computational time compared to a conventional 2D FE analysis. This reduction in computational time comes at the cost of an error that, for a certain level of accuracy, cannot be reduced further by refining the 1D mesh, but can be reduced by computing higher order semi-analytic approximations (see Appendix A). The semi-analytic approximation was initially derived for an infinite domain, but edge effects can be included. This is particularly useful for the efficient simulation of repetitive scanning patterns that are typically encountered in lithography and additive manufacturing applications. In the presented wafer heating example, the first-order semi-analytic approximation reduces the time to compute the temperature field resulting from the scanning of a single field by a factor 10 compared to a standard FE approach with similar 4%-accuracy. For the second-order semi-analytic approximation this error is reduced to 1.5%. However, a much larger amount of CPU time is saved by assembling the response for the real expose pattern using superposition, see Fig. 12, where the conventional FE approach leads to excessive CPU times.
As stated before, the proposed method introduces an additional error that at a certain accuracy level cannot be decreased by refining the 1D mesh size. However, the 4% error achieved for the firstorder semi-analytic approximation in the wafer heating example is acceptable for two reasons: (1) the temperature field itself is not of interest, but only the displacements induced by it. For the considered example, the 4% error in the first-order semi-analytic approximation of the temperature field leads to a 0.4% error in the resulting displacements; (2) efficient simulation of the moving heat load problem is especially important when the model is used in a feedback control loop or for model-based controller design. In these cases, the feedback controller creates some robustness to modeling errors.
For the considered example, the difference between the 4%accuracy of the first-order semi-analytic approximation and the 1.5%-accuracy of the second-order semi-analytic approximation is relatively small. So, for the considered rectangular uniform shape of the heat load considering more terms in the Taylor series approximation (11) does not rapidly increase the accuracy of the semi-analytic approximation. On the other hand, for an element size of L e ¼ 0:8 mm, the second-order approximation is still computed 10 times faster than an FE solution with similar accuracy, so that higher-order approximations may be valuable. Moreover, for a smoother shape XðxÞ of the heat load in the x-direction increasing the order of the semi-analytic approximation may lead to faster convergence.
In the derivation of the approximation, we assumed that the heat load is of the form (2), with XðxÞ P 0; YðyÞ P 0, and Q ðtÞ P 0. The assumption that XðxÞ; YðyÞ, and Q ðtÞ are nonnegative is not very restrictive. Consider for example the case where Q ðtÞ does not satisfy this assumption. In that case we can always write Q ðtÞ as the difference of two nonnegative functions QðtÞ ¼ Q þ ðtÞ À Q À ðtÞ; ð44Þ where Q þ ðtÞ ¼ Q ðtÞ when Q ðtÞ > 0 and zero otherwise, and Q À ðtÞ ¼ À Q ðtÞ when Q ðtÞ < 0 and zero otherwise. Since the functions Q þ ðtÞ and Q À ðtÞ are nonnegative, we can compute the approximation when we replace Q ðtÞ in (1) by Q þ ðtÞ and by Q À ðtÞ. Since the PDE in (3) is linear, we may substract the response for Q À ðtÞ from the response for Q þ ðtÞ to find the response for Q ðtÞ. A similar procedure can be applied when XðxÞ or YðyÞ are not nonnegative.
Finally, we note that the method can be generalized to a 3D spatial domain ðx; y; zÞ 2 R 3 . It is shown in Appendix B that the problem can be decoupled into four 1D problems, when the applied heat load can be written as Q 3D ðx; y; z; tÞ ¼ XðxÞYðy À vtÞZðzÞ Q ðtÞ: ð45Þ  In this case, the potential reduction in computational cost is even larger than the reduction for the two-dimensional problem considered in this paper. Uðz 0 ; sÞZðz À z 0 Þ dz 0 :