A New Approach to Solve Non-Fourier Heat Equation via Empirical Methods Combined with the Integral Transform Technique in Finite Domains

Abstract

This chapter deals with the validity/limits of the integral transform technique on finite domains. The integral transform technique based upon eigenvalues and eigenfunctions can serve as an appropriate tool for solving the Fourier heat equation, in the case of both laser and electron beam processing. The crux of the method consists in the fact that the solutions by mentioned technique demonstrate strong convergence after the 10 eigenvalues iterations, only. Nevertheless, the method meets with difficulties to extend to the case of non-Fourier equations. A solution is however possible, but it is bulky with a weak convergence and requires the use of extra-boundary conditions. To surpass this difficulty, a new mix approach is proposed with this chapter resorting to experimental data, in order to support a more appropriate solution. The proposed method opens in our opinion a beneficial prospective for either laser or electron beam processing.

1. Introduction

1.1 Mathematical background

The heat equation can be solved in a simpler mode via the Fourier heat equation, which involves the propagation of heat waves with infinite speed. This hypothesis is in particular valid for many applications, such as laser-metal interaction in the frame of two-temperature model [1, 2].

The solution of Fourier equations can be inferred using different mathematical techniques via Green function, integral, Laplace transform, or complex analysis. The predictions of the solutions given by the mentioned methods are of analytical or semianalytical nature and confirm the experimental data for certain situations such as laser–metal interaction.

One basically assumes that the heat waves propagation speed is inversely proportional to the square root of the relaxation time. A smaller relaxation time leads to higher heat speed waves, resulting in a good Fourier approximation. If one requires however a more accurate description of experimental data, one should introduce a more exact method to solve the non-Fourier equation involving a finite heat wave speed.

A mixed solution of the non-Fourier equation combines the theoretical method of finite integral transforms with information from experimental data. Thus, two additional boundary conditions can be imposed, which will lead to a semianalytical solution of the non-Fourier equation. The finite domains of the integral transform method for Fourier equations are eigenfunctions and values, which reach after 10 iterations a quite conform solution for the Fourier equation [3, 4, 5, 6]. This method is applied to the non-Fourier equation, and the final form is obtained, with the support of experimental results.

A new heat transfer model was adopted in order to unify the thermal field distribution in both laser and electron beam processing. An analytical solution using non-Fourier heat equation has been developed corresponding to boundary conditions in the case of material processing. The model has been compared with the experimental data obtained using an in-house developed facility. A simplified and easy-to-use model via MATHEMATICA software stands for the novelty of the current work.

2. Non-Fourier equation

The non-Fourier equation is hyperbolic and can be written as:

$$\frac{{\partial }^{2}T}{\partial r^2}+\frac{1}{r}\frac{\partial T}{\partial r}+\frac{{\partial }^{2}T}{\partial z^2}+\frac{1}{r^2}\frac{{\partial }^{2}T}{\partial {\phi }^2}-\frac{1}{\gamma }\frac{\partial T}{\partial t}-\frac{{\tau }_0}{\gamma }\frac{{\partial }^{2}T}{\partial t^2}=-\frac{P\left(r,\phi ,z,t\right)}{K}.$$E1

Here, T is target temperature, r is the radial coordinate, z is the spatial coordinate on the direction of laser beam propagation, t is time, ${\mathit{\tau }}_{\mathbf{0}}$ is relaxation time, and $\mathit{\gamma }$ is thermal diffusivity. P stands for the source term,$\mathit{\phi }$ is the angular coordinate, while K is the target thermal conductivity. For a simple general solution, one assumes a cylindrical target with angular symmetry, under irradiation with a Gaussian laser beam with the center at r = z = 0. In this case, the temperature does not depend on $\mathit{\phi }$:

$$\frac{\partial T}{\partial \phi }=0.=>T=T\left(r,z,t\right).$$E2

The corresponding boundary conditions are:

$$K\frac{\partial T\left(r,z,t\right)}{\partial r}{│}_{r=b}+\mathit{hT}\left(b,z,t\right)=0.$$E3

Here r = b is the cylinder radius, while h is the heat transfer coefficient. The boundary conditions for the z coordinate are:

$$K\frac{\partial T\left(r,z,t\right)}{\partial z}{│}_{z=0}-\mathit{hT}\left(r,0,t\right)=0,$$E4

and

$$K\frac{\partial T\left(r,z,t\right)}{\partial z}{│}_{z=a}+\mathit{hT}\left(r,a,t\right)=0,$$E5

where a is the cylinder length. We will pass on the effective solution of the equation. The operator D_r was defined as:

$$D_r=\frac{{\partial }^2}{\partial r^2}+\frac{1}{r}\left(\frac{\partial }{\partial r}\right).$$E6

This applies to Eq. (1) with the boundary conditions:

$$\frac{{\partial }^2{K}_r}{\partial r^2}+\frac{1}{r}\frac{\partial K_r}{\partial r}+{\mu }^2{K}_r=0,$$E7

and

$$K\frac{\partial K_r\left(r,z,t\right)}{\partial r}{│}_{r=b}+h{K}_r\left(b,z,t\right)=0.$$E8

Eqs. (7) and (8) corroborate to Eq. (9):

$$\frac{h}{K}{J}_0\left({\mu }_{i}b\right)-{\mu }_i{J}_1\left({\mu }_{i}b\right)=0.$$E9

One can deduce, based upon the theory of finite integral transforms, the eigenfunction ${\stackrel{\sim }{K}}_r\left(r,{\mu }_i\right)$ corresponding to the eigenvalue ${\mu }_i$:

$${\stackrel{\sim }{K}}_r\left(r,{\mu }_i\right)=J_O\left({\mu }_{i}r\right)r\frac{1}{C_i},$$E10

where the normalization constant is given by:

$$C_i\left(r\right)={\int }_0^{b}r{K}_r^2\left(r,{\mu }_i\right)\mathit{dr}=\frac{b^2}{2{\mu }_i^2}\left(\frac{h^2}{k^2}+{\mu }_i^2\right)J_0^2\left({\mu }_{i}b\right).$$E11

One defines:

$$\stackrel{\sim }{T}\left({\mu }_{i,}z,t\right)=\frac{1}{C_i}{\int }_0^{b}T\left(r,z,t\right)r{J}_0\left({\mu }_{i}r\right)\mathit{dr},$$E12

and

$$\stackrel{\sim }{P}\left({\mu }_i,z,t\right)=\frac{1}{C_i}{\int }_0^{b}P\left(r,z,t\right)r{J}_0\left({\mu }_{i}r\right)\mathit{dr}.$$E13

Eq. (1) in this case converts to:

$$-{\mu }_i^2·\stackrel{\sim }{T}+\frac{{\partial }^2\stackrel{\sim }{T}}{\partial z^2}-\frac{1}{\gamma }\frac{\partial \stackrel{\sim }{T}}{\partial t}-\frac{{\tau }_0}{\gamma }\frac{{\partial }^2\stackrel{\sim }{T}}{\partial t^2}=-\frac{\stackrel{\sim }{P}}{K}.$$E14

One obtains for z coordinate via similar mathematical calculation:

$$\frac{{\partial }^2{K}_z}{\partial z^2}+{\lambda }^2{K}_z=0,$$E15

and

$${\left[K\frac{\partial K_Z}{\partial z}-h{K}_z\right]}_{z=0}=0,$$E16

as well as:

$${\left[K\frac{\partial K_Z}{\partial z}+h{K}_z\right]}_{z=a}=0.$$E17

The − and + signs for h in Eqs. (16) and (17) denote the target heat absorption and emission, respectively. One has:

$$K_z\left(z,\lambda \right)=\mathit{cos}\left(\lambda z\right)+\frac{h}{\mathit{\lambda K}}\mathit{sin}\left(\lambda z\right),$$E18

and

$$2cot\left({\lambda }_{j}a\right)=\frac{{\lambda }_{j}k}{h}-\frac{h}{{\lambda }_{j}k}.$$E19

Here, ${\mathit{\lambda }}_{\mathit{j}}$ denotes eigenvalues along the z-axis. From theory [7], it follows:

$${\stackrel{\sim }{K}}_z\left(z,{\lambda }_j\right)=\frac{1}{C_j}{K}_z\left(z,{\lambda }_j\right),$$E20

and

$$C_j={\int }_0^a{K}_z^2\left(z,{\lambda }_j\right)\mathit{dz}.$$E21

Note that Eqs. (12) and (13) discuss the eigenvalues along the r-axis, only. After introducing the eigenvalues along the z-axis, one can step ahead to generalize Eqs. (12) and (13) as:

$$\overline{T}\left({\mu }_i,{\lambda }_{j,}t\right)=\frac{1}{C_i{C}_j}{\int }_0^a{\int }_0^{b}T\left(r,z,t\right)r{K}_r\left({\mu }_i,r\right)K_z\left({\lambda }_j,z\right)\mathit{\text{drdz}},$$E22

and

$$\overline{P}\left({\mu }_i,{\lambda }_j,t\right)=\frac{1}{C_i{C}_j}{\int }_0^a{\int }_0^{b}P\left(r,z,t\right)r{K}_r\left({\mu }_i,r\right)K_z\left({\lambda }_j,z\right)\mathit{\text{drdz}}.$$E23

Eq. (1) becomes now:

$$\left({\mu }_i^2+{\lambda }_j^2\right)\overline{T}+\frac{1}{\gamma }\frac{\partial \overline{T}}{\partial t}+\frac{{\tau }_0}{\gamma }\frac{{\partial }^2\overline{T}}{\partial t^2}=\frac{\overline{P}}{K}.$$E24

We next applied the direct and inverse Laplace integral transform to solve Eq. (1) in relation to time. C[1] and C[2] stand for the normalizing coefficients with respect to the experimental data. The results are as follows:

$$T\left(r,z,t,{\tau }_0\right)=\sum _{i=1}^{10}\sum _{j=1}^{10}\left[\begin{array}{c}\frac{P\left({\mu }_i,{\lambda }_j\right)}{{\mu }_i^2+{\lambda }_j^2}+C\left[1\right]{℮}^{\frac{\left(-\frac{1}{\gamma }-\sqrt{\frac{1}{{\gamma }^2}-4{\mu }_i^2\frac{{\tau }_0}{\gamma }-4{\lambda }_j^2\frac{{\tau }_0}{\gamma }}\right)\bullet t}{\frac{2{\tau }_0}{\gamma }}}+C\left[2\right]\\ {℮}^{\frac{\left(-\frac{1}{\gamma }+\sqrt{\frac{1}{{\gamma }^2}-4{\mu }_i^2\frac{{\tau }_0}{\gamma }-4{\lambda }_j^2\frac{{\tau }_0}{\gamma }}\right)t}{\frac{2{\tau }_0}{\gamma }}}\end{array}\right]\left(K_r\left({\mu }_i,r\right)K_z\left({\lambda }_{j)},z\right)\right).$$E25

and

$$T\left({\mu }_i,{\lambda }_j\right)=\frac{1}{C_i{C}_j}{\int }_0^a{\int }_0^{b}P\left(r,z,t\right)r{K}_r\left({\mu }_i,r\right)K_z\left({\lambda }_j,z\right)\mathit{\text{drdz}}.$$E26

We finally mention that for an intermediate point in the experimental curve, one has:

$$T\left(r,z,t\right)=T\left({\tau }_0,C\left[1\right],C\left[2\right],r,z,t\right).$$E27

With the boundary conditions:

$$T\left(r,z,0,C\left[1\right],C\left[2\right],{\tau }_0\right)={22}^{°}\mathrm{C}.$$E28

$$T\left(r,z,\infty ,C\left[1\right],C\left[2\right],{\tau }_0\right)={22}^{°}\mathrm{C}.$$E29

3. Two-temperature model in the non-Fourier version

The two-temperature model (TTM) is based upon two coupled equations:

$$A{T}_e\left(\frac{\partial T_e}{\partial t}\right)+\frac{K{\tau }_0}{\gamma }\left(\frac{{\partial }^2{T}_e}{\partial t^2}\right)=K\left(\frac{{\partial }^2{T}_e}{\partial x^2}+\frac{{\partial }^2{T}_e}{\partial y^2}+\frac{{\partial }^2{T}_e}{\partial z^2}\right)-G\left(T_e-T_i\right)+P_a\left(\overrightarrow{r},t\right),$$E30

$$C_i\left(\frac{\partial T_i}{\partial t}\right)=G\left(T_e-T_i\right).$$E31

Here $T_e$ and $T_i$ stand for the electron and phonon temperatures, respectively. G is the coupling factor between electrons and phonons. $P_a\left(\overrightarrow{r},t\right)$ is the heat source, which is induced via laser-metal interaction. The interaction could be considered of either classical or steady-state quantum mechanical type. A is the electron heat capacity, and K is the thermal conductivity of the metal. According to Ref. [8], $G$ can be determined from:

$$G=\frac{{\pi }^2\mathit{mN}{v}^2}{6\tau T_i}{\left(\frac{T_e}{T_i}\right)}^4\times {\int }_0^{\raisebox{1ex}{$T_e$}\!\left/ \!\raisebox{-1ex}{$T_d$}\right.}\left[x^4/\left(e^x-1\right)\right]\mathit{dx},$$E32

where m is the electron mass, N is the conduction electron density, v is the velocity of sound in the metal, $\tau$ is the electron–phonon collision time, and $T_D$ is the Debye temperature. The exact data for each metal (Cu, Ag, Al, or Fe) are available from text books and current literature (e.g., [8, 9, 10]). As for all metals, $C_i>>A$, one may assume, in a first approximation, that:

$$K\left(\frac{{\partial }^2{T}_e}{\partial x^2}+\frac{{\partial }^2{T}_e}{\partial y^2}+\frac{{\partial }^2{T}_e}{\partial z^2}\right)-C_i\frac{\partial T_i}{\partial t}-\frac{K{\tau }_0}{\gamma }\left(\frac{{\partial }^2{T}_e}{\partial t^2}\right)=-P_a\left(\overrightarrow{r},t\right).$$E33

According to the Nolte model [2], one has:

$$T_i=\kappa T_e,$$E34

where

$$\kappa =\frac{{\tau }_L}{{\tau }_L+{\tau }_i}.$$E35

Here, ${\tau }_i$ is the lattice cooling time while ${\tau }_L$ is the pulse duration time. Consequently, one has:

$$\frac{\partial T_i}{\partial t}=\kappa \frac{\partial T_e}{\partial t}.$$E36

Eq. (3) can be rewritten as:

$$K\left(\frac{{\partial }^2{T}_e}{\partial x^2}+\frac{{\partial }^2{T}_e}{\partial y^2}+\frac{{\partial }^2{T}_e}{\partial z^2}\right)-C_i\kappa \frac{\partial T_e}{\partial t}-\frac{K{\tau }_0}{\gamma }\left(\frac{{\partial }^2{T}_e}{\partial t^2}\right)=-P_a\left(\overrightarrow{r},t\right).$$E37

It follows that:

$$\left(\frac{{\partial }^2{T}_e}{\partial x^2}+\frac{{\partial }^2{T}_e}{\partial y^2}+\frac{{\partial }^2{T}_e}{\partial z^2}\right)-\frac{1}{\gamma }\frac{\partial T_e}{\partial t}-\frac{{\tau }_0}{\gamma }\left(\frac{{\partial }^2{T}_e}{\partial t^2}\right)=-\frac{-P_a\left(\overrightarrow{r},t\right)}{K}$$E38

with

$$\gamma =\frac{K}{C_i\kappa }.$$E39

Under the most general form, the heat source reads as:

$$P_a={\sum }_{m,n}{I}_{\mathit{mn}}\left(y,z\right)\left({\alpha }_{\mathit{mn}}{e}^{-{\alpha }_{\mathit{mn}}x}\right)\left(1-r_{\mathit{Smn}}\right)+r_{\mathit{Smn}}\delta \left(x\right)+q\left.c\right)\left(H\left(t\right)-H\left(t-t_0\right)\right).$$E40

Here, $I_{\mathit{mn}}\left(y,z\right)$ stands for the laser transverse mode {m,n} while ${\alpha }_{\mathit{mn}}$ is the linear absorption coefficient, and $r_{\mathit{Smn}}$ is the surface absorption coefficient. The quantum corrections (qc) are steady state for the respective mode. H stands for the step function, and t₀ for the exposure time. One model explains the continuous laser beam irradiation, while the other one, more realistic in our opinion, illustrates the laser beam in pulse form. The equivalence between the two models requires therefore that the intensity versus time plot should cover the same area.

In order to make a comparison with experiments, one needs besides analytical description, concrete numerical values. The next step is therefore to estimate the eigenvalues, numerically. For this purpose, Eq. (7) can be solved using the integral transform technique, and eigenfunctions and eigenvalues could be calculated. One has three differential equations as follows (K represents the eigenfunctions, while $\lambda ,\mu ,\text{and}\xi$ are the eigenvalues) [5]:

$$\frac{{\partial }^2{K}_x}{\partial x^2}+{\lambda }_i^2{K}_x=0,$$E41

$$\frac{{\partial }^2{K}_y}{\partial y^2}+{\mu }_j^2{K}_y=0,$$E42

$$\frac{{\partial }^2{K}_z}{\partial z^2}+{\xi }_k^2{K}_z=0.$$E43

The final solutions could be achieved on the basis of Eqs. (41)–(43):

$$K_x=\mathit{cos}\left({\lambda }_{i}x\right)+\frac{h}{K{\lambda }_i}\mathit{sin}\left({\lambda }_{i}x\right),$$E44

$$K_y=\mathit{cos}\left({\mu }_{j}y\right)+\frac{h}{K{\mu }_j}\mathit{sin}\left({\mu }_{j}y\right),$$E45

$$K_z=\mathit{cos}\left({\xi }_{k}z\right)+\frac{h}{K{\xi }_k}\mathit{sin}\left({\xi }_{k}z\right).$$E46

The boundary conditions are:

$${\left[\frac{\partial K_x}{\partial x}-\frac{h{K}_x}{K}\right]}_{x=0}=0;{\left[\frac{\partial K_x}{\partial x}+\frac{h{K}_x}{K}\right]}_{x=a}=0,$$E47

$${\left[\frac{\partial K_y}{\partial y}+\frac{h{K}_y}{K}\right]}_{y=0}=0;{\left[\frac{\partial K_y}{\partial y}+\frac{h{K}_y}{K}\right]}_{y=b}=0,$$E48

$${\left[\frac{\partial K_z}{{\partial }_z}+\frac{h{K}_z}{K}\right]}_{z=0}=0;{\left[\frac{\partial K_z}{\partial z}+\frac{h{K}_z}{K}\right]}_{z=c}=0.$$E49

Here, a, b, and c are the metal sample dimensions. As for the boundary conditions, the eigenvalues (h is the heat transfer coefficient) can be inferred from Eqs. (47) to (49), as:

$$2\mathit{cot}\left({\lambda }_{i}a\right)=\frac{{\lambda }_{i}K}{h}-\frac{h}{K{\lambda }_i},$$E50

$$2\mathit{cot}\left({\mu }_{j}b\right)=\frac{{\mu }_{j}K}{h}-\frac{h}{K{\mu }_j},$$E51

$$2\mathit{cot}\left({\xi }_{k}c\right)=\frac{{\xi }_{k}K}{h}-\frac{h}{K{\xi }_k}.$$E52

The solution is obtained via integral transform technique as:

$$\begin{array}{l}T\left(x,y,z,t,C\left[1\right],C\left[2\right],{\tau }_0\right)\\ =\sum _{i=1}^{10}\sum _{j=1}^{10}\sum _{k=1}^{10}\left[\begin{array}{c}\frac{P\left({\mu }_i,{\lambda }_{j,}{\xi }_k\right)}{{\mu }_i^2+{\lambda }_j^2+{\xi }_k^2}+C\left[1\right]{℮}^{\frac{\left(-\frac{1}{\gamma }-\sqrt{\frac{1}{{\gamma }^2}-4{\mu }_i^2\frac{{\tau }_0}{\gamma }-4{\lambda }_j^2\frac{{\tau }_0}{\gamma }-4{\xi }_k^2\frac{{\tau }_0}{\gamma }}\right)\bullet t}{\frac{2{\tau }_0}{\gamma }}}+C\left[2\right]\\ {℮}^{\frac{\left(-\frac{1}{\gamma }+\sqrt{\frac{1}{{\gamma }^2}-4{\mu }_i^2\frac{{\tau }_0}{\gamma }-4{\lambda }_j^2\frac{{\tau }_0}{\gamma }-4{\xi }_k^2\frac{{\tau }_0}{\gamma }}\right)\bullet t}{\frac{2{\tau }_0}{\gamma }}}\end{array}\right]\\ \times \left(K_x\left({\mu }_i,x\right)K_y\left({\lambda }_j,y\right)K_z\left({\xi }_k,z\right)\right)\end{array}$$E53

The advantage of Eq. (53) in our model is related to a quick converging series. Thus, after 10 iterations, the solution’s accuracy reaches already 10⁻² K in the case of thermal distribution [11].

4. Experimental details

The experimental setup is operated by a Nd:YAG pulsed laser source (λ = 355 nm) (Surelite II from Continuum), generating pulses of 6 ns duration with (130 ± 0.6) mJ energy at a frequency repetition rate of 10 Hz. The laser beam had a spatial top-hat distribution. The laser beam was focused onto the metallic target surface by a lens with 240 mm focal lengths. An Al bulk target of (10 × 10 × 5) mm³ was used in experiments. The laser fluence was set at ∼7.5 J/cm² to surpass the ablation threshold but also to avoid the excessive plasma formation. A crater of 18 μm depth was dig into the sample after the application of 1000 subsequent laser pulses, as checked up by a Vernier Caliper instrument. During the multipulse laser irradiation, the thermal distribution was monitored on the sample back side via a FLUKA thermocouple connected to a computer having Lab view software, while the sample was irradiated at the top. All experiments were performed on an in-house developed equipment at Laser Department, National Institute for Laser, Plasma and Radiation Physics (INFLPR), Magurele, Romania.

5. Results and discussion

Experiments and simulations were carried out during the heating of a metallic target. The boundary conditions were described by Eq. (28). In all figures, the experimental data are plotted with dots while the simulations are represented by a continuous line. Relaxation time,${\tau }_0,$ was assumed 0.5 ps (Figure 1), 1 ns (Figure 2), and 1 μs (Figure 3), respectively. For simulation, a heat transfer coefficient = 3 × 10⁻⁷ W mm⁻² K⁻¹ was selected. As known [12, 13, 14], for a very low heat transfer coefficient, the eigenvalues are positive very small numbers, resulting in a linear thermal distribution curve, as visible in Figures 1–3.

Figure 1.

Time evolution of temperature for a relaxation time of 0.5 ps: experiments (dotted line) vs. simulation (continuous line).

Figure 2.

Time evolution of the temperature for a relaxation time of 1 ns: experiments (dotted line) vs. simulation (continuous line).

Figure 3.

Time evolution of the temperature for a relaxation time of 1 μs: experiments (dotted line) vs. simulation (continuous line).

The best agreement between theory and experiment was achieved for a relaxation time,${\tau }_0,$of 0.5 ps, as visible from Figure 1. We note that this is in accordance with available literature on subject [15].

6. Conclusions and outlook

The two-temperature model was generalized to the case of the non-Fourier approach via the electron-phonon relaxation time. Boundary conditions, Eq. (28) for heating and Eq. (29) for cooling, were considered to this purpose. The obtained solutions prove useful for experimental data analysis. The mathematical method belongs to the eigenvalues and functions family, while details on software are available from Ref. [6].

The exact nature of the metallic target (in our case aluminum) could be detected from the electron-phonon relaxation time using integral transform technique mix via acquired experimental data. The method can be extended to any experimental sample (metal) with the high accuracy.

Acknowledgments

CNM, MO, NM, and CR acknowledge for financial support by Romanian Ministry of Research, Innovation and Digitalization, under Romanian National NUCLEU Program LAPLAS VI–contract no. 16N/2019. CNM, NM, CR, and INM thank for the financial support from a grant of the Romanian Ministry of Education and Research, CNCS-UEFISCDI, project number ID code RO-NO-2019-0498 and UEFISCDI under the TE_196/2021 and PED_306/2020. M.A.M. received financial support from the European Union’s Horizon 2020 (H2020) research and innovation program under the Marie Skłodowska–Curie grant agreement no. 764935.