Sage Journals: Discover world-class research

Abstract

The Rosenthal steady-state analytical solution for the temperature distribution caused by a moving point heat source on a semi-infinite, homogeneous, isotropic solid has been extensively used in modeling metallurgical processes, e.g., arc welding. This study develops a three-dimensional analytical closed-form solution for the temperature field induced by a point heat source moving across a semi-infinite orthotropic solid. The formulation accommodates arbitrary orientations of the heat source’s motion relative to the material’s in-plane principal axes and is extended to solids with finite thickness. Subsequently, using the superposition of the linear solutions, a general methodology is proposed to predict the temperature distribution resulting from an arbitrarily distributed heat source. Verification of the closed-form solution and validation of the distributed heating condition against finite element simulations demonstrate excellent agreement. The analytical framework offers potential for thermal modeling of processing methods used for polymeric composites, e.g., additive manufacturing and continuous welding of thermoplastic composites.

Keywords

Analytical solution moving heat source orthotropic solid three-dimensional thermal analysis

Introduction

With the widespread integration of automation in the industrial processing of various materials and the manufacturing of structural components, the theoretical study of moving heat sources in localized processing has garnered attention since the 1940s with the advent of automated metal welding, and has remained highly relevant with recent advancements in additive manufacturing processes for metallic alloys. Understanding heat propagation within a material subjected to a moving heat source is critical for predicting its thermal history, which directly governs the evolution of microstructural properties and impacts the final product’s quality in terms of mechanical performance and dimensional stability.

Rosenthal¹ developed a foundational theory for thermal analysis of a moving point heat source on isotropic solids. He provided analytical quasi-steady solutions for one-dimensional, two-dimensional, and three-dimensional heat flow in infinite and semi-infinite solids bounded by planes. Meanwhile, Jaeger² and Carslaw and Jaeger³ made seminal contributions in the development of the analytical solutions of the heat conduction equation for different moving source shapes (e.g., band, square, and rectangular with uniform distribution). Later, Hou and Komanduri⁴ extended the solutions for additional shapes (e.g., elliptical and circular geometries) and other heat distributions (e.g., parabolic and normal profiles). Despite the advancements, the original Jaeger-Rosenthal closed-form analytical solution for a moving point heat source has remained the well-established framework in the process modelling of joining or additive manufacturing (e.g., beam welding, arc welding, selective laser melting, and laser metal deposition) of isotropic metallic materials (e.g., stainless steel, aluminum alloys, titanium alloys),⁵ expressing the temperature distribution in a semi-infinite solid as:

T = T_{\infty} + \frac{Q^{'}}{2 π k r} e^{- \frac{U}{2 γ} (r + x)},

(1)

where

T_{\infty}

is the temperature at an infinite distance from the source,

Q^{'}

is the heat power transferred to the semi-infinite solid,

k

is thermal conductivity,

γ

is thermal diffusivity,

U

is magnitude of the source velocity,

r

is the absolute distance from the source, and

x

is the distance from the source projected along the velocity vector, as depicted in Figure 1.

Figure 1.

Representation of the thermal field via moving reference frame attached to a moving point heat source on semi-infinite isotropic solid⁵.

While analytical modeling of a moving heat source on isotropic materials is well-established and no longer considered novel, there remains a need for closed-form analytical solutions for arbitrarily oriented moving heat sources on orthotropic materials, where the thermal conductivity is orientation-dependent. Developing a theoretical framework for orthotropic materials would assist in modeling the thermal processing of fiber-reinforced composites, with examples of polymer-based composites listed in Figure 2. In automated fiber placement (AFP) or automated tape laying (ATL) of thermoset composites (TSCs), typically, an infrared heating system integrated into the placement head increases the temperature of the substrate for enhanced tackiness as the head moves forward during deposition. In AFP of thermoplastic composites (TPCs), where diode laser,⁶ hot gas torch,⁷ or flash lamp⁸ heating systems are typically used to heat the material above its melting point for in-situ consolidation, or for fused filament fabrication (FFF) of fiber reinforced polymers (FRPs) as a similar example, the temperature distribution around the deposition area controls the quality of the final part. A lack of precise prediction over temporal and spatial temperature distribution within the material can lead to violating the process window and causing excessive defects such as material degradation, incomplete healing, and increased void content,⁹ which ultimately compromise the mechanical performance of the structure.^10,11 Robust and rapid predictive thermal tools allow for online fine-tuning of process parameters, such as the speed, intensity, or distribution of the moving heat source, to achieve the optimal material properties while minimizing energy usage and maximizing efficiency.

Figure 2.

Application of moving heat source in different processes of polymeric composites.

Among the original studies on thermal modeling for in-situ consolidation of thermoplastic composites, Tierney and Gillespie¹² employed a one-dimensional steady-state model to predict the through-thickness temperature distribution when the composite surface is exposed to hot gas impingement. Ghasemi Nejhad et al.¹³ developed a two-dimensional analytical model for the thermal analysis of thermoplastic tape laying, modeling a moving uniform heat source over a transversely isotropic domain. Pitchumani et al.¹⁴ also utilized a two-dimensional thermal model to predict the evolution of interfacial bonding, polymer degradation, and void consolidation, with the aim of defining the optimal processing window. While one- and two-dimensional analytical and numerical models continue to be efficient and popular,¹⁵ their accuracy might be limited – for instance, by neglecting heat diffusion across the width. This has motivated the adoption of three-dimensional numerical approaches, especially finite element methods,¹⁶ enabling more accurate thermal predictions, though at the cost of substantially increased computational effort.

The requirement for accurate thermal models is equally critical in continuous in-situ joining techniques for thermoplastic composites (e.g., induction welding,¹⁷ ultrasonic welding,¹⁸ friction stir welding¹⁹) where localized melting and solidification of interfacial thermoplastics occurs within the order of a second during the forward motion of the end-effector. On the other hand, in cutting or surface treatment of composites,²⁰ predicting temperature distribution caused by direct heat source (e.g., laser beam, thermal spray, respectively) or friction-induced heat generation (e.g., CNC milling, automated sanding, respectively) may allow for prevention of thermally softened or minimization of degraded material by optimizing process parameters.

The application of an analytical framework for orthotropic solids can even be extended to additively manufactured isotropic materials if significant porosity, characterized by longitudinal voids/gaps along the deposition axis, renders the homogenized properties effectively orthotropic. Beyond polymers and composites, the framework may be relevant to specific metallic substrates exhibiting orthotropic thermal properties. An example could be single-crystal alloys, such as those used in turbine blades,²¹ which display direction-dependent thermal conduction and are subjected to laser-based repair techniques involving localized, moving heat sources.

In this study, an analytical closed-form solution is developed for the temperature distribution induced by a point heat source moving over a semi-infinite orthotropic solid, considering an arbitrary orientation of motion relative to the material’s principal axes. The solution is then extended to accommodate solids with finite thickness. Finally, leveraging the principle of superposition for linear systems, a methodology is proposed to predict the temperature distribution resulting from an arbitrarily distributed heat source. The proposed solutions for the point heat source (in both semi-infinite and finite-thickness cases) and the results for the distributed heat source are validated against finite element simulations, demonstrating excellent agreement.

Theoretical modelling

The following equation represents the three-dimensional heat conduction for an orthotropic solid in Cartesian coordinates:

k_{x} \frac{\partial^{2} T}{\partial {x^{'}}^{2}} + k_{y} \frac{\partial^{2} T}{\partial {y^{'}}^{2}} + k_{z} \frac{\partial^{2} T}{\partial {z^{'}}^{2}} + q = ρ c_{p} \frac{\partial T}{\partial t} .

(2)

The spatial coordinates $x^{'}$ - $y^{'}$ - $z^{'}$ are defined with respect to a stationary reference frame fixed to the material and aligned along its principal axes. Temperature ( $T$ ) can be a function of spatial coordinates and time, while the principal thermal conductivities ( $k_{x}, k_{y}, k_{z}$ ), specific heat capacity ( $c_{p}$ ), and density ( $ρ$ ) are considered constant properties of the homogeneous domain. The term $q$ represents the rate of internal heat generation per unit volume. In the case of an infinitesimal point heat source, $q$ is nonzero only at the precise location of the heat source and is equal to zero elsewhere throughout the domain. Hence, for a single point heat source located at ( $x_{0}^{'}, y_{0}^{'}$ , $z_{0}^{'}$ ) at time $t$ , $q$ can be represented as:

q = Q (t) δ (x^{'} - x_{0}^{'}) δ (y^{'} - y_{0}^{'}) δ (z^{'} - z_{0}^{'}),

(3)

where

Q

denotes instantaneous thermal power injected into the system by the point heat source. For now, since the interest is to analyze the heat conduction behavior in the domain excluding the singular point of the heat source, it is evident that

q = 0

Consider a point heat source moving horizontally ( $v_{z} = 0$ ) inside an infinite solid domain at a constant velocity $v$ . It will be shown that this particular choice reflects the configuration of the most relevant applications – such as additive manufacturing and continuous welding – where the end-effector traverses a horizontal path. A moving reference frame is defined such that it is fixed to the point heat source but remains aligned with the principal axes of the orthotropic material, resulting in a coordinate transformation to the $x$ - $y$ - $z$ system, as illustrated in Figure 3.

Figure 3.

Representation of stationary $x^{'}$ - $y^{'}$ - $z^{'}$ reference frame, and moving $x$ - $y$ - $z$ reference frame attached to the moving point heat source inside the infinite orthotropic solid (x-axis oriented along the principal axis of the material, e.g., fibers).

Let $θ$ denote the angle between the velocity vector v and the positive x-axis. Accordingly, the components of the velocity vector are $v_{x} = v \cos (θ)$ and $v_{y} = v \sin (θ)$ . Now, the relationship between the stationary and moving coordinate systems is described by:

x = x^{'} - v_{x} t,

(4)

y = y^{'} - v_{y} t,

(5)

z = z^{'} .

(6)

To mathematically simplify the orthotropic nature, spatial coordinates and velocity components can be transformed into scaled variables so that they are normalized based on the corresponding thermal conductivities:

X = \frac{x}{\sqrt{k_{x}}},

(7)

Y = \frac{y}{\sqrt{k_{y}}},

(8)

Z = \frac{z}{\sqrt{k_{z}}},

(9)

X^{'} = \frac{x^{'}}{\sqrt{k_{x}}},

(10)

Y^{'} = \frac{y^{'}}{\sqrt{k_{y}}},

(11)

Z^{'} = \frac{z^{'}}{\sqrt{k_{z}}},

(12)

V_{x} = \frac{v_{x}}{\sqrt{k_{x}}},

(13)

V_{y} = \frac{v_{y}}{\sqrt{k_{y}}} .

(14)

With these, equations (4), (5), and (6) transform to:

X = X^{'} - V_{x} t,

(15)

Y = Y^{'} - V_{y} t,

(16)

Z = Z^{'} .

(17)

By the transformation of partial derivatives using the chain rule, the quasi-steady heat equation can be expressed in the new coordinates as follows:

\frac{\partial^{2} T}{\partial X^{2}} + \frac{\partial^{2} T}{\partial Y^{2}} + \frac{\partial^{2} T}{\partial Z^{2}} = - ρ c_{p} V_{x} \frac{\partial T}{\partial X} - ρ c_{p} V_{y} \frac{\partial T}{\partial Y} .

(18)

Inspired by Rosenthal,¹ to mathematically separate the effects of advection and symmetric heat diffusion, and to further simplify equation (18), consider:

T - T_{\infty} = e^{- \frac{1}{2} ρ c_{p} V_{x} X} e^{- \frac{1}{2} ρ c_{p} V_{y} Y} φ (X, Y, Z),

(19)

where

T_{\infty}

represents the temperature of the solid far from the influence of the heat source, where temperature gradients are effectively zero. The first and second terms on the right-hand side describe the exponential decay in the positive X and Y directions, respectively, as a direct result of the heat source moving forward. The function

φ (X, Y, Z)

, which accounts for the effective symmetric heat diffusion, must be determined to satisfy equation (18). By substituting equation (19) into equation (18) and applying the necessary manipulations (detailed derivations are shown in Appendix A):

\frac{\partial^{2} φ}{\partial X^{2}} + \frac{\partial^{2} φ}{\partial Y^{2}} + \frac{\partial^{2} φ}{\partial Z^{2}} - {(\frac{1}{2} ρ c_{p} \sqrt{{V_{x}}^{2} + {V_{y}}^{2}})}^{2} φ = 0 .

(20)

This equation appears symmetric with respect to the coordinates, and technically, $φ$ is invariant of orientation. In this case, let us define a new parameter:

s = \sqrt{X^{2} + Y^{2} + Z^{2}},

(21)

where

s

accounts for the radius of a sphere in the

X

Y

Z

coordinates (or an ellipsoid in the

x

y

z

coordinates). Meanwhile, the Laplacian in spherical coordinates (

s

α

β

) appears as:

\nabla^{2} = \frac{1}{s^{2}} \frac{\partial}{\partial s} (s^{2} \frac{\partial}{\partial s}) + \frac{1}{s^{2} \sin α} \frac{\partial}{\partial α} (\sin α \frac{\partial}{\partial α}) + \frac{1}{s^{2} \sin^{2} α} \frac{\partial^{2}}{\partial β^{2}} .

(22)

Since $φ$ is radially symmetric ( $φ = φ (s)$ ), the derivatives with respect to $α$ and $β$ vanish. Therefore, by transforming the Laplacian from Cartesian to spherical coordinates, the equation (20) transforms into:

\frac{1}{s^{2}} \frac{\partial}{\partial s} (s^{2} \frac{\partial φ}{\partial s}) - {(\frac{1}{2} ρ c_{p} \sqrt{{V_{x}}^{2} + {V_{y}}^{2}})}^{2} φ = 0 .

(23)

Let us define

g (s) = s φ (s)

. Then, equation (23) becomes:

\frac{\partial^{2} g}{\partial s^{2}} - {(\frac{1}{2} ρ c_{p} \sqrt{{V_{x}}^{2} + {V_{y}}^{2}})}^{2} g = 0 .

(24)

The solution of this linear second-order homogeneous differential equation is:

g = C_{1} e^{- \frac{1}{2} ρ c_{p} \sqrt{{V_{x}}^{2} + {V_{y}}^{2}} s} + C_{2} e^{+ \frac{1}{2} ρ c_{p} \sqrt{{V_{x}}^{2} + {V_{y}}^{2}} s},

(25)

φ = \frac{C_{1}}{s} e^{- \frac{1}{2} ρ c_{p} \sqrt{{V_{x}}^{2} + {V_{y}}^{2}} s} + \frac{C_{2}}{s} e^{+ \frac{1}{2} ρ c_{p} \sqrt{{V_{x}}^{2} + {V_{y}}^{2}} s} .

(26)

However, the physics of the problem requires

C_{2}

to be zero, since the influence of the heat source must decay to zero at infinity. By applying this and substituting equation (26) into equation (19):

T = T_{\infty} + e^{- \frac{1}{2} ρ c_{p} (V_{x} X + V_{y} Y)} \frac{C_{1}}{s} e^{- \frac{1}{2} ρ c_{p} \sqrt{{V_{x}}^{2} + {V_{y}}^{2}} s} .

(27)

Up to this point, the boundary condition at infinity has been applied, leading to $C_{2} = 0$ . To determine $C_{1}$ , it is necessary to apply the remaining boundary condition, which is the singular point of the heat source where the heat power $Q$ is injected. This requires incorporating the term $q$ back into the analysis, as the singular point of the heat source is now of primary interest. By applying a transformation to equation (3), converting from the $x^{'}$ - $y^{'}$ - $z^{'}$ coordinates to the $X$ - $Y$ - $Z$ coordinates:

q = \frac{Q}{\sqrt{k_{x}} \sqrt{k_{y}} \sqrt{k_{z}}} δ (X) δ (Y) δ (Z) .

(28)

To simplify the determination of the constant $C_{1}$ , a special case is considered where the velocity of the point heat source is zero. Under this condition, firstly, equation (27) simplifies to:

T = T_{\infty} + \frac{C_{1}}{s} .

(29)

Secondly, the heat equation in the transformed coordinates may appear as:

\frac{\partial^{2} T}{\partial X^{2}} + \frac{\partial^{2} T}{\partial Y^{2}} + \frac{\partial^{2} T}{\partial Z^{2}} = - \frac{Q}{\sqrt{k_{x}} \sqrt{k_{y}} \sqrt{k_{z}}} δ (X) δ (Y) δ (Z) .

(30)

Let us define $P = \frac{Q}{\sqrt{k_{x}} \sqrt{k_{y}} \sqrt{k_{z}}}$ , assigning it artificial dimension of [W]. With this definition, equation (30) resembles heat propagation from a point source with heat power $P$ in an isotropic solid, assuming a thermal conductivity of $k = 1$ [W/m⋅K] in the $X$ - $Y$ - $Z$ coordinate system. Also, based on Fourier’s law of heat conduction in the radial direction:

{\bar{q}}_{s} = - k \frac{\partial T}{\partial s} .

(31)

Here, ${\bar{q}}_{s}$ denotes the rate of heat conduction per unit area in the radial direction. In the particular scenario of heat propagation from a point source within an infinite isotropic solid (i.e., purely spherical diffusion), the radial heat flux corresponds to the input heat power divided by the area of a spherical surface centered at the point of heat injection. Based on this, and considering equation (29), equation (31) can be reformulated as:

\frac{P}{4 π s^{2}} = \frac{C_{1}}{s^{2}} .

(32)

Therefore:

C_{1} = \frac{Q}{4 π \sqrt{k_{x}} \sqrt{k_{y}} \sqrt{k_{z}}} .

(33)

With the constant

C_{1}

determined, the general form of equation (27) becomes:

T = T_{\infty} + \frac{Q}{4 π \sqrt{k_{x} k_{y} k_{z}} s} e^{- \frac{1}{2} ρ c_{p} (V_{x} X + V_{y} Y + \sqrt{{V_{x}}^{2} + {V_{y}}^{2}} s)} .

(34)

By substituting equations (7), (8), (9), (13), (14), and (21) into equation (34), the temperature distribution in the $x$ - $y$ - $z$ coordinate system can be expressed as:

T = T_{\infty} + \frac{Q}{4 π \sqrt{k_{x} k_{y} k_{z}} \sqrt{\frac{x^{2}}{k_{x}} + \frac{y^{2}}{k_{y}} + \frac{z^{2}}{k_{z}}}} e^{- \frac{1}{2} ρ c_{p} v (\frac{x \cos θ}{k_{x}} + \frac{y \sin θ}{k_{y}} + \sqrt{\frac{\cos^{2} θ}{k_{x}} + \frac{\sin^{2} θ}{k_{y}}} \sqrt{\frac{x^{2}}{k_{x}} + \frac{y^{2}}{k_{y}} + \frac{z^{2}}{k_{z}}})} .

(35)

It is often more appropriate to represent the temperature distribution in the $\bar{x}$ - $\bar{y}$ - $\bar{z}$ coordinate system, which shares the same origin as the $x$ - $y$ - $z$ system, and $\bar{z}$ remains coincident with $z$ , while $\bar{x}$ is aligned along the velocity vector. In this scenario, a two-dimensional rotational coordinate transformation is applied to equation (35):

T = T_{\infty} + \frac{Q}{4 π \sqrt{k_{x} k_{y} k_{z}} \sqrt{\frac{{(\bar{x} \cos θ - \bar{y} \sin θ)}^{2}}{k_{x}} + \frac{{(\bar{x} \sin θ + \bar{y} \cos θ)}^{2}}{k_{y}} + \frac{{\bar{z}}^{2}}{k_{z}}}} \times e^{- \frac{1}{2} ρ c_{p} v (\frac{(\bar{x} \cos θ - \bar{y} \sin θ) \cos θ}{k_{x}} + \frac{(\bar{x} \sin θ + \bar{y} \cos θ) \sin θ}{k_{y}} + \sqrt{\frac{\cos^{2} θ}{k_{x}} + \frac{\sin^{2} θ}{k_{y}}} \sqrt{\frac{{(\bar{x} \cos θ - \bar{y} \sin θ)}^{2}}{k_{x}} + \frac{{(\bar{x} \sin θ + \bar{y} \cos θ)}^{2}}{k_{y}} + \frac{{\bar{z}}^{2}}{k_{z}}})} .

(36)

The presented solution provides insight into the temperature distribution within an infinite solid, in which a moving point source with power $Q$ travels, and where temperature is solely a function of spatial coordinates. This steady-state condition holds since the reference frame is attached to the heat source moving with constant velocity. The temperature at any given point remains constant in this moving frame, where the point is stationary with respect to the frame but moves relative to the material. However, when considering the thermal history of the material – where time-dependent properties (e.g., cooling-rate-dependent crystallization of thermoplastic polymers²²) become relevant – the instantaneous rate of temperature change on the material can still be expressed within the moving reference frame, as detailed in Appendix B:

\dot{T} = \frac{T - T_{\infty}}{\frac{{(\bar{x} \cos θ - \bar{y} \sin θ)}^{2}}{k_{x}} + \frac{{(\bar{x} \sin θ + \bar{y} \cos θ)}^{2}}{k_{y}} + \frac{{\bar{z}}^{2}}{k_{z}}} \times [v (\frac{(\bar{x} \cos θ - \bar{y} \sin θ) \cos θ}{k_{x}} + \frac{(\bar{x} \sin θ + \bar{y} \cos θ) \sin θ}{k_{y}}) + \frac{1}{2} ρ c_{p} v^{2} \times (\frac{(\bar{x} \cos θ - \bar{y} \sin θ) \cos θ}{k_{x}} + \frac{(\bar{x} \sin θ + \bar{y} \cos θ) \sin θ}{k_{y}}) \times \sqrt{\frac{\cos^{2} θ}{k_{x}} + \frac{\sin^{2} θ}{k_{y}}} \times \sqrt{\frac{{(\bar{x} \cos θ - \bar{y} \sin θ)}^{2}}{k_{x}} + \frac{{(\bar{x} \sin θ + \bar{y} \cos θ)}^{2}}{k_{y}} + \frac{{\bar{z}}^{2}}{k_{z}}} + \frac{1}{2} ρ c_{p} v^{2} \times (\frac{\cos^{2} θ}{k_{x}} + \frac{\sin^{2} θ}{k_{y}}) \times (\frac{{(\bar{x} \cos θ - \bar{y} \sin θ)}^{2}}{k_{x}} + \frac{{(\bar{x} \sin θ + \bar{y} \cos θ)}^{2}}{k_{y}} + \frac{{\bar{z}}^{2}}{k_{z}})] .

(37)

The surface centerline ( $\bar{y} = 0, \bar{z} = 0$ ) is often of particular interest in practice for evaluating representative temperature distribution and representative heating and cooling rates. This condition simplifies equations (36) and (37) to the following forms, respectively:

T = T_{\infty} + \frac{Q}{4 π \sqrt{k_{x} k_{y} k_{z}} \sqrt{\frac{\cos^{2} θ}{k_{x}} + \frac{\sin^{2} θ}{k_{y}}} \bar{x}} e^{- ρ c_{p} v (\frac{\cos^{2} θ}{k_{x}} + \frac{\sin^{2} θ}{k_{y}}) \bar{x}},

(38)

\dot{T} = (T - T_{\infty}) \times [\frac{v}{\bar{x}} + ρ c_{p} v^{2} (\frac{\cos^{2} θ}{k_{x}} + \frac{\sin^{2} θ}{k_{y}})] .

(39)

Returning to equation (36), while the developed closed-form solution represents the temperature distribution resulting from a point heat source moving within an infinite orthotropic solid, due to the symmetry of heat propagation with respect to the $\bar{x}$ - $\bar{y}$ plane ( $\bar{z} = 0$ ), the same formulation can be applied to model the temperature distribution induced by a point heat source injecting power $Q^{'} = Q / 2$ , moving along the adiabatic surface of a semi-infinite orthotropic solid. Therefore, the term $Q$ can be simply replaced by $2 Q^{'}$ to adapt the solution for a moving heat source with power $Q^{'}$ on the surface of a semi-infinite solid. The verification of this analytical model using finite element analysis will be subsequently discussed.

Orthotropic Solids with finite Thickness

The next step involves evaluating the effect of a moving heat source on solids with finite thickness. For constant and known velocity, material properties, and orientation, as shown in equation (36), the temperature distribution can be represented by the elementary function $f$ :

T - T_{\infty} = f (Q, \bar{x}, \bar{y}, \bar{z}) .

(40)

As a linear system, the temperature distribution in the solid domain induced by multiple heat sources can be obtained by summation of the effects of all individual sources. Consider an infinite number of point heat sources, spaced $2 h$ apart along the $\bar{z}$ axis, each with identical power and velocity vectors, moving through an infinite solid as depicted in Figure 4. The resulting temperature distribution is given by:

T - T_{\infty} = \sum_{κ = - K}^{K} f (Q, \bar{x}, \bar{y}, \bar{z} + 2 κ h), K \to \infty .

(41)

Figure 4.

Vertically aligned multiple identical point heat sources and corresponding planes of symmetry.

In this configuration, two sets of symmetry planes emerge, functioning as adiabatic boundaries due to the zero heat flux across symmetry planes. The first set consists of horizontal planes passing through the point heat source, one of which is depicted by a dashed line in Figure 4. The second set includes horizontal planes midway between two adjacent heat sources, one of which is shown as a dotted line in Figure 4. The domain $- h \leq \bar{z} \leq 0$ can effectively represent a solid of finite thickness with adiabatic top and bottom surfaces, on which a moving point heat source with power $Q / 2$ injects heat into the system at the location ( $\bar{x} = 0, \bar{y} = 0, \bar{z} = 0$ ).

Two-dimensional Solution

Although the primary focus of this study is not on two-dimensional planar heat transfer, such analysis becomes relevant in certain applications – such as thin plates – where through-thickness thermal gradients are insignificant and in-plane heat conduction dominates. An exact two-dimensional solution cannot be obtained as a limiting case of equation (36). Instead, the governing equations should be derived independently, as presented in Appendix C. The final expression for the temperature field is:

T = T_{\infty} + \frac{\bar{Q}}{2 π \sqrt{k_{x} k_{y}}} [e^{- \frac{1}{2} ρ c_{p} v (\frac{(\bar{x} \cos θ - \bar{y} \sin θ) \cos θ}{k_{x}} + \frac{(\bar{x} \sin θ + \bar{y} \cos θ) \sin θ}{k_{y}})}] \times K_{0} (\frac{1}{2} ρ c_{p} v \sqrt{\frac{\cos^{2} θ}{k_{x}} + \frac{\sin^{2} θ}{k_{y}}} \sqrt{\frac{{(\bar{x} \cos θ - \bar{y} \sin θ)}^{2}}{k_{x}} + \frac{{(\bar{x} \sin θ + \bar{y} \cos θ)}^{2}}{k_{y}}}) .

(42)

Here,

\bar{Q}

denotes the heat source power per unit depth [W/m], and

K_{0}

represents the modified Bessel function of the second kind of order zero.

Verification

Finite element analysis was performed using COMSOL Multiphysics, employing quadratic Lagrange elements for temperature discretization. In the first step, to perform verification of the analytical model for a semi-infinite solid, the model geometry consisted of a rectangular solid block with dimensions 0.2 m × 0.05 m × 0.01 m. These dimensions were chosen large enough to ensure that the region of interest around the point heat source remained unaffected by boundary effects, thereby satisfying the assumption of a semi-infinite solid. Boundary conditions were specified as follows: constant room temperature was applied to the two lateral faces, the bottom face, and the front (upstream) face, while thermal insulation was imposed on the top face and the back (downstream) face. A point heat source was positioned at the center of the top face (origin) to inject a constant heat power into the system. The relative motion of the heat source was incorporated by defining a translational motion for the solid domain with velocity field of −v in the $\bar{x}$ direction to enable steady-state analysis. In this configuration, the point heat source was treated as stationary within the computational domain, while capturing a relative velocity of v in the $\bar{x}$ direction with respect to the solid.

The mesh was constructed using tetrahedral elements. Mesh parameters included a maximum element size of 2 × 10⁻³ m, a minimum element size of 1 × 10⁻⁶ m, and a maximum growth rate of 1.05. This meshing strategy resulted in a total of 1,373,879 elements, with the highest mesh resolution concentrated in the vicinity of the point heat source to capture critical thermal gradients, while coarser resolution was employed in regions distant from the area of interest to optimize computational efficiency (see Figure 5 (top)). The input parameters of the simulation are as follows, where the selected material as our example is Carbon-Fiber/Poly-Ether-Ether-Ketone (CF/PEEK) (Table 1):

Figure 5.

Finite element model in COMSOL Multiphysics.

Table 1.

Values of CF/PEEK thermal properties²³ and parameters used for FEM verification.

Parameters	Values
$k_{x}$	5.9 W/(m*K)
$k_{y}$	0.7 W/(m*K)
$k_{z}$	0.7 W/(m*K)
$ρ$	1575 kg/(m^3)
$c_{p}$	1300 J/(kg*K)
$v$	1 cm/s
$Q^{'}$	2 W
$T_{\infty}$	20°C
$θ$	0°, 45°, 90°

Figure 5 (middle) and (bottom) present color maps of the temperature distribution calculated using the FEM in the $\bar{x}$ - $\bar{y}$ and $\bar{x}$ - $\bar{z}$ planes, respectively. These results confirm that the block’s width and thickness were sufficiently large to replicate heat propagation in a semi-infinite solid within the region of interest around the point heat source.

The results of the analytical model for a semi-infinite solid and the finite element analysis are compared in Figure 6, represented as isothermal contours for three different orientations of the velocity vector relative to the material’s principal x-axis. For each orientation, temperature distributions in the $\bar{x}$ - $\bar{y}$ and $\bar{x}$ - $\bar{z}$ planes are presented. The FEM results are shown as discrete data points, while the analytical solution is depicted using continuous colored lines. The excellent agreement between the two approaches confirms the validity of equation (36). The contours indicate that when the velocity vector aligns with one of the material’s principal axes (e.g., aligned with x ( $θ = 0 °$ ) or with y ( $θ = 90 °$ )), the temperature distribution exhibits symmetry with respect to the $\bar{x}$ - $\bar{z}$ plane. In contrast, when the heat source is moving off-axis (e.g., θ = 45°), the temperature distribution becomes asymmetric relative to the $\bar{x}$ - $\bar{z}$ plane.

Figure 6.

(a). Temperature distributions from analytical solution and FEM for a semi-infinite solid with heat source velocity oriented at 0° to the material x-axis: in-plane (top) and through-thickness (bottom). (b). Temperature distributions from analytical solution and FEM for a semi-infinite solid with heat source velocity oriented at 45° to the material x-axis: in-plane (top) and through-thickness (bottom). (c). Temperature distributions from analytical solution and FEM for a semi-infinite solid with heat source velocity oriented at 90° to the material x-axis: in-plane (top) and through-thickness (bottom).

In another stage of verification, the analytical model for a solid with finite thickness (equation (41)) was compared with finite element analysis (FEA). For this case, the geometry of the FEA model was modified by significantly reducing the thickness to 0.25 × 10⁻³ m to account for the influence of the bottom face on the temperature distribution. The boundary condition at the bottom face was then set to adiabatic, while all other parameters, including those listed in Table 1, remained unchanged. The comparison of the analytical and FEA results for the in-plane and through-thickness temperature distributions are illustrated in Figure 7 for three different orientations of the material’s x-axis relative to the heat source velocity. Once again, the excellent agreement between the analytical and FEA results confirms the validity of the analytical framework. In the analytical calculations for this case, the parameter K in equation (41) was set to only 5 (superposition of 11 heat sources), which proved sufficient to achieve high accuracy in the region of interest around the heat source. Due to the exponential decay of the heat source’s effect with distance, selecting a large value for K is unnecessary to maintain acceptable accuracy.

Figure 7.

(a). Temperature distributions from analytical solution and FEM for a finite-thickness solid with heat source velocity oriented at 0° to the material x-axis: in-plane (top) and through-thickness (bottom). (b). Temperature distributions from analytical solution and FEM for a finite-thickness solid with heat source velocity oriented at 45° to the material x-axis: in-plane (top) and through-thickness (bottom). (c). Temperature distributions from analytical solution and FEM for a finite-thickness solid with heat source velocity oriented at 90° to the material x-axis: in-plane (top) and through-thickness (bottom).

Distributed heat source

The theoretical framework for thermal analysis using a moving point heat source is only valid in practical applications when the heating system delivers highly localized power to the surface. This imposes a limitation on the direct implementation of the existing analytical solution for practical systems. To overcome this constraint, the surface area of the distributed heat flux, which effectively governs the temperature distribution within the solid, is discretized into a series of point heat sources. This discretization aims to replicate the thermal influence of a continuous heat flux distribution. Leveraging the linearity of the thermal system, the temperature at any point within the solid domain can be computed as the superposition of the temperature increments (ΔT) induced by each individual point heat source.

Let $p (\bar{x}, \bar{y})$ represent the given continuous heat flux profile, with the origin of $\bar{x}$ - $\bar{y}$ - $\bar{z}$ coordinate system positioned at the peak of $p$ . The effective surface area of the heat flux can be discretized into a grid of $(M_{f} + M_{b} + 1) \times (N_{l} + N_{r} + 1)$ point heat sources, spaced at intervals of $Δ ξ$ in the $\bar{x}$ -direction and $Δ η$ in the $\bar{y}$ -direction, as illustrated in Figure 8.

Figure 8.

Representation of a continuous heat source as an array of discrete point sources.

The temperature at an arbitrary excited point (ep) located at $(\bar{x}, \bar{y}$ , $\bar{z})$ within an orthotropic solid with an adiabatic bottom surface can be determined by summing the temperature increments induced by all source points (sp). This can be expressed as:

T - T_{\infty} = \sum_{m = - M_{b}}^{M_{f}} \sum_{n = - N_{r}}^{N_{l}} \sum_{κ = - K}^{K} f (p_{(m ∆ ξ, n ∆ η)} Δ ξ Δ η, \bar{x} - m ∆ ξ, \bar{y} - n ∆ η, \bar{z} + 2 κ h) .

(43)

The bounds of the summation indices, $m$ (ranging from $- M_{b}$ to $M_{f}$ ) and $n$ (ranging from $- N_{r}$ to $N_{l}$ ), which define the computational domain, need to be determined with an understanding of the heat flux distribution profile. The spatial extent of the summation, represented as a rectangular area with dimensions $(M_{f} + M_{b}) Δ ξ \times (N_{l} + N_{r}) Δ η$ , should be sufficiently large to capture the region responsible for at least 99% of the total heat exchange with the surface.

An additional consideration in defining the computational domain is the choice of $Δ ξ$ and $Δ η$ , which need to be sufficiently small to avoid significant truncation errors where steep gradients of the heat distribution exist. Another critical consideration in selecting the sizes of $Δ ξ$ and $Δ η$ arises when it comes to the location of target excited point for temperature calculation. When the excited point is located closer to the surface, the discrete source points must be positioned closer together (i.e., smaller $Δ ξ$ and $Δ η$ ) to more effectively replicate the effects of a continuous heating profile. This requirement will be clarified in the following example of calculations.

Even though the heat flux profile $p (\bar{x}, \bar{y})$ could be any arbitrary function, including non-elementary functions (e.g., a piecewise function), or discrete sets of data points, in this example we adopt a 2D Gaussian function with elliptical decay as the heat source which can approximate a laser beam or a light source:

p = p_{\max} e^{- (\frac{{\bar{x}}^{2}}{a^{2}} + \frac{{\bar{y}}^{2}}{b^{2}})} .

(44)

The constants

p_{\max}

a

, and

b

selected in the upcoming analysis are listed in Table 2. The thickness and material orientation are also presented in Table 2, while the selected material properties and the velocity of the heat source are the same as those presented in Table 1.

Table 2.

Parameters and corresponding values used for FEM verification.

Parameters	Values
$p_{\max}$	6E5 W/(m^ 2)
$a$	5E-3 m
$b$	2E-3 m
$h$	2.5E-3 m
$θ$	90°

Calculations were performed at 3000 points within the solid, along a line segment spanning from (10 mm, 0 mm, −0.2 mm) to (−30 mm, 0 mm, −0.2 mm), for various heat source spacings ( $Δ ξ$ = $Δ η$ ) of 0.1 mm, 0.2 mm, 0.4 mm, 0.8 mm, 1.2 mm, 1.6 mm, and 2.0 mm. The resulting temperature distributions, as shown in Figure 9, indicate that when the source spacing becomes too large relative to the depth of the target point ( $\bar{z})$ , significant fluctuations emerge. These fluctuations, which stem from the nature of discrete point sources, deviate from the continuous nature of the actual heat flux profile. An analysis across different source spacings, depths, and heat flux intensities (Figure 9 presents results for a single depth only) suggests that for achieving high accuracy, defined as an amplitude of fluctuation waves less than 0.1% of local magnitude of ( $T - T_{\infty})$ , the following condition must be satisfied for the spacing size:

\max (\frac{Δ ξ}{| \bar{z} |}, \frac{Δ η}{| \bar{z} |}) \leq 1 .

(45)

Figure 9.

Predicted temperature distribution along the axis of motion for different point source spacing ratios.

In Figure 10, the temperature distributions predicted by the proposed model are shown at various depths within the $\bar{x} - \bar{z}$ plane and $\bar{y} - \bar{z}$ plane and are compared with finite element simulation results with a continuous heat flux profile on the surface using equation (44) and parameters in Table 2. To maintain computational efficiency, the spatial discretization ratios are set as $\frac{Δ ξ}{| \bar{z} |} = \frac{Δ η}{| \bar{z} |} = 1$ . It is worthy of attention that direct computation of the temperature distribution at the surface ( $\bar{z} = 0$ ) is not feasible since $Δ ξ$ and $Δ η$ must approach zero to achieve sufficient resolution. To overcome this, the surface temperature can be determined by extrapolating temperature values from deeper layers. For example, linear extrapolation using data at $\bar{z} = - d$ and $\bar{z} = - d / 2$ can approximate the temperature at $\bar{z} = 0$ . The accuracy of the extrapolation depends on the choice of $d$ , which needs to be selected small enough to ensure capturing the variations effectively for high accuracy at regions with steep gradient. A practical guideline for the size of $d$ is to ensure that:

\frac{| T (\bar{x}, \bar{y}, - \frac{d}{2}) - T (\bar{x}, \bar{y}, - d) |}{T (\bar{x}, \bar{y}, - \frac{3 d}{4}) - T_{\infty}} < 10 % .

(46)

Figure 10.

Temperature distribution along (left plot) and perpendicular (right plot) to the direction of source motion: Validation of the model for a distributed heat source against FEM results at different depths.

For example, the predicted surface temperature results shown in Figure 10 are obtained by extrapolating the temperature curves at $\bar{z} = - 0.1$ mm and $\bar{z} = - 0.05$ mm to $\bar{z} = 0$ . In the presented results, the error in the predicted surface temperature using extrapolation is below 2%, whereas the directly computed temperatures for deeper layers exhibit errors below 0.1%.

The steady-state simulation computation time for the COMSOL finite element solver on a typical personal computer was measured at 38 seconds, while the computation time for each point using the analytical model was on the order of centi-seconds (e.g., 0.029 seconds for $\frac{Δ ξ}{| \bar{z} |}$ = 1) on the same computer. This can highlight the utility of the analytical approach for rapid temperature estimations when the purpose is to determine temperature at a specific point or along lines consisting of a moderate number of points. The finite element solver in COMSOL solves a large system of algebraic equations derived from the discretized heat transfer equations using advanced iterative methods, such as the GCRO-DR (Generalized Conjugate Residual with Deflated Restarts) method, being specifically optimized for computational efficiency in solving large linear systems, making it well-suited for the current linear simulation. However, while the finite element method requires the definition of boundary conditions over a domain that extends far beyond the region of interest to ensure accurate results, the analytical approach evades this dependency. The analytical solution directly computes the temperature at desired locations without being influenced by the solution for adjacent points or the requirement of defining far-field conditions. The independence from boundary condition definitions in remote regions and the ability to compute temperature directly at a target point, simplifies the analytical approach and enhances its practicality. In the meantime, the small computation time may show potential for real-time thermal predictions when it comes to smart process control systems in continuous composite manufacturing/joining processes, e.g., AFP or induction welding, which enables dynamic adjustments to process parameters to maintain operations within the optimal window.

A parametric analysis comparing different fiber orientations and heat source velocities under the same distributed heating condition (Table 2) and material properties (Table 1) is presented in Figure 11. The results are visualized using the iso-surface corresponding to the critical glass transition temperature (Tg = 143 $℃$ ) of CF/PEEK, with the maximum width, length, and depth of the heat-soaked region annotated for each case.

Figure 11.

Iso-surface of T = 143 $℃$ for different heat source velocities and material orientations depicted with relative [ $\bar{x}$ :1, $\bar{y}$ :1, $\bar{z}$ :10] scale of magnification (yellow point marking the origin of the reference frame, i.e., the center of the heat source profile).

The results indicate that at a velocity of 1 cm/s, reducing the fiber angle from 90° to 0° leads to a moderate reduction (8%) in the width of the region exceeding Tg , while the length of this region increases substantially (80%). A similar trend is observed at 2 cm/s, where the same variation in fiber orientation results in a 5% reduction in width and a 54% increase in length. Regarding the maximum penetration depth of the heated region, decreasing the fiber angle from 90° to 0° results in a 55% increase at 1 cm/s and a 43% increase at 2 cm/s. These findings suggest that as the axis of maximum thermal conductivity (i.e., typically fiber direction for CFRPs) becomes more parallel to the direction of heat source motion, the thermally affected zone becomes not only longer but also deeper. Furthermore, at lower velocities, the dimensions of the heated zone exhibit greater sensitivity to fiber orientation than at higher velocities.

From a different standpoint, increasing the heat source velocity from 1 cm/s to 2 cm/s reduces the size of the region, as the material is exposed to the heat source for a shorter duration. The reduction in length is quantified as 13%, 17%, and 26% for fiber angles of 90°, 45°, and 0°, respectively, while the corresponding reductions in width are 20%, 19%, and 17%. The most pronounced effect is observed in depth reduction, which decreases by 43%, 44%, and 48% for fiber angles of 90°, 45°, and 0°, respectively. Overall, while this analysis highlights the influence of fiber orientation and velocity on the dimensions of the heat-soaked region, the sensitivity of these effects is directly proportional to the material’s degree of orthotropy (e.g., $k_{x} / k_{y}$ ) in case of a general material.

Conclusion

This paper presents a closed-form solution for the thermal analysis of a moving point heat source over a semi-infinite orthotropic solid with arbitrary orientation (e.g., off-axis motion). The resulting temperature distribution is expressed in equation (36), while the local rate of temperature change in the material is formulated in equation (37). While equation (36) was derived for a point heat source of power $Q$ moving within an infinite solid, it remains directly applicable to the case of a moving heat source of power $Q^{'} = Q / 2$ on the surface of a semi-infinite solid, as a result of symmetry considerations. The solution is then extended to account for finite-thickness solids with adiabatic boundaries, broadening its applicability to practical scenarios. The accuracy of the proposed solution has been verified through very close agreement with results from finite element analysis. Furthermore, a discretization methodology is introduced to approximate continuous heat sources with arbitrary shapes and distributions by superimposing multiple point heat sources, which has also been validated against the finite element method to prove its excellent reliability. The proposed theoretical framework offers a computationally efficient approach for rapid thermal predictions, making it applicable to simulating thermal distributions in the processing of orthotropic composite materials. Finally, although the primary focus of this work was the development of a three-dimensional analytical framework, a two-dimensional solution was also derived and succinctly presented in equation (42), offering potential applicability to thin-plate scenarios.

Footnotes

Acknowledgements

Financial support from Natural Sciences and Engineering Research Council of Canada is appreciated.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: Natural Sciences and Engineering Research Council of Canada.

ORCID iD

Mahmoud Fereidouni

Appendix

References

Rosenthal

. The theory of moving sources of heat and its application to metal treatments. Trans Am Soc Mech Eng 1946; 68(8): 849–865.

Jaeger

. Moving sources of heat and the temperature at sliding contacts. Proc. Roy. Soc. New South Wales 1942; 76: 203.

Carslaw

Jaeger

. Conduction of heat in solids. 2nd ed. Oxford: Oxford University Press, 1959.

Hou

Komanduri

. General solutions for stationary/moving plane heat source problems in manufacturing and tribology. Int J Heat Mass Tran 2000; 43(10): 1679–1698.

Moda

Chiocca

Macoretta

, et al. Technological implications of the Rosenthal solution for a moving point heat source in steady state on a semi-infinite solid. Mater Des 2022; 223: 110991.

Hudson

Dolph

Grose

, et al. Thermal response of thermoplastic composite tape during in-situ consolidation automated fiber placement using a laser heat source. SAMPE Conference and Exhibition, 2023.

Tafreshi

Hoa

Shadmehri

, et al. Heat transfer analysis of automated fiber placement of thermoplastic composites using a hot gas torch. Adv Manuf Polym Compos Sci 2019; 5(4): 206–223.

Tavkari

Ganesan

Amgai

, et al. Thermal characterization of Xenon-Arc flash lamp heating system for automated fiber placement (AFP) of thermoplastic composites. Compos Appl Sci Manuf 2025; 190: 108568.

Fereidouni

Hoa

. In-situ consolidation of thermoplastic composites by automated fiber placement: characterization of defects. J Thermoplast Compos Mater 2025; 38(2): 471–511.

10.

Yap

Heathman

Bidabadi

, et al. In-plane properties of an in-situ consolidated automated fiber placement thermoplastic composite. Compos Appl Sci Manuf 2025; 188: 108525.

11.

Rastak

Van Hoa

. Effect of fiber orientation on the strength of thermoplastic composite tubes subject to four-point bending. Composites Part C: Open Access 2025; 16: 100569.

12.

Tierney

Gillespie

JJW

. Modeling of heat transfer and void dynamics for the thermoplastic composite tow-placement process. J Compos Mater 2003; 37(19): 1745–1768.

13.

Ghasemi Nejhad

Cope

Güceri

. Thermal analysis of in-situ thermoplastic composite tape laying. J Thermoplast Compos Mater 1991; 4(1): 20–45.

14.

Pitchumani

Gillespie

JJW

Lamontia

. Design and optimization of a thermoplastic tow-placement process with in-situ consolidation. J Compos Mater 1997; 31(3): 244–275.

15.

Stelter

Sreekantamurthy

Hudson

, et al. Thermal modelling of the in-situ consolidation of automated fiber placement of thermoplastic composites. InSociety for the Advancement of Material and Process Engineering (SAMPE) 2024.

16.

Zacherl

Fontes

Shadmehri

. Three-dimensional heat transfer analysis of hot gas torch (HGT)-assisted automated fiber placement (AFP) for thermoplastic composites. Compos Struct 2024; 343: 118256.

17.

Lionetto

Zecca

Buccoliero

, et al. Experimental and numerical investigation of susceptor-aided continuous induction welding of low-melt PAEK composites. Polym Compos 2025.

18.

Larsen

Endrass

Jarka

, et al. Exploring ultrasonic and resistance welding for thermoplastic composite structures: process development and application potential. Compos B Eng 2025; 289: 111927.

19.

Kumar

Medhi

Roy

. Friction stir welding of thermoplastic composites. In Advances in Industrial and Production Engineering: Select Proceedings of FLAME. Singapore: Springer Singapore, 2019, vol 24, pp. 221–235.

20.

Chahine

Marathe

Collins

, et al. Effect of plasma treatment on LMPAEK/CF tape and composites manufactured by automated tape placement (ATP). Compos Appl Sci Manuf 2025; 188: 108540.

21.

Swanson

Arakere

. Effect of crystal orientation on analysis of single-crystal, nickel-based turbine blade superalloys NASA report. NASA 2000: 65, TP-2000-210074.

22.

Hoa

. Principles of the manufacturing of composite materials. 2nd ed., DEStech Publications, 2018.

23.

Cogswell

. Thermoplastic aromatic polymer composites: a study of the structure, processing and properties of carbon fibre reinforced polyetheretherketone and related materials. Elsevier, 2013.

24.

Abramowitz

Stegun

. Handbook of mathematical functions: with formulas, graphs, and mathematical tables. Courier Corporation 1965.

An analytical approach for three-dimensional thermal analysis of moving heat source on orthotropic solid

Abstract

Keywords

Introduction

Theoretical modelling

Orthotropic Solids with finite Thickness

Two-dimensional Solution

Verification

Distributed heat source

Conclusion

Footnotes

Acknowledgements

Declaration of conflicting interests

Funding

ORCID iD

Appendix

References