Next Article in Journal
Birnbaum-Saunders Quantile Regression Models with Application to Spatial Data
Next Article in Special Issue
Lp-Solution to the Random Linear Delay Differential Equation with a Stochastic Forcing Term
Previous Article in Journal
Modeling an Uncertain Productivity Learning Process Using an Interval Fuzzy Methodology
Previous Article in Special Issue
Two New Strategies for Pricing Freight Options by Means of a Valuation PDE and by Functional Bounds
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Second-Order Dual Phase Lag Equation. Modeling of Melting and Resolidification of Thin Metal Film Subjected to a Laser Pulse

1
Faculty of Mechanical Engineering, Silesian University of Technology, 44-100 Gliwice, Poland
2
Department of Technical Sciences, University of Occupational Safety Management, 40-007 Katowice, Poland
*
Author to whom correspondence should be addressed.
Mathematics 2020, 8(6), 999; https://doi.org/10.3390/math8060999
Submission received: 2 June 2020 / Revised: 15 June 2020 / Accepted: 16 June 2020 / Published: 18 June 2020
(This article belongs to the Special Issue Models of Delay Differential Equations)

Abstract

:
The process of partial melting and resolidification of a thin metal film subjected to a high-power laser beam is considered. The mathematical model of the process is based on the second-order dual phase lag equation (DPLE). Until now, this equation has not been used for the modeling of phase changes associated with heating and cooling of thin metal films and the considerations regarding this issue are the most important part of the article. In the basic energy equation, the internal heat sources associated with the laser action and the evolution of phase change latent heat are taken into account. Thermal processes in the domain of pure metal (chromium) are analyzed and it is assumed that the evolution of latent heat occurs at a certain interval of temperature to which the solidification point was conventionally extended. This approach allows one to introduce the continuous function corresponding to the volumetric fraction of solid or liquid state at the neighborhood of the point considered, which significantly simplifies the phase changes modeling. At the stage of numerical computations, the authorial program based on the implicit scheme of the finite difference method (FDM) was used. In the final part of the paper, the examples of numerical computations (including the results of simulations for different laser intensities and different characteristic times of laser pulse) are presented and the conclusions are formulated.

1. Introduction

Heat transfer through thin films subjected to an ultrafast laser pulse is of vital importance in microtechnology applications and is a reason that the problems related to the fast heating of solids have become a very active research area. The problems of melting/resolidification processes modeling, which may be the result of heating with a laser beam, are also important from the technical point of view. So far, the method using the equation based on the second order model with two delay times for the phase changes modeling was not presented in literature. This was the most important motivation for the authors to undertake research in this area. In Section 5, the comparison of the results obtained with the similar solution on the basis of the first-order dual phase lag equation (DPLE) is presented.
The mathematical model of macroscale heat conduction is based on the parabolic Fourier equation. This equation was formulated under the assumption of instantaneous propagation of the thermal wave in the domain considered. It is obvious that this assumption is not correct, but for the problems concerning the analysis of macroscale heat conduction processes, the obtained results are fully satisfactory. Despite this, the attempts have been made to modify the Fourier equation to a form that better reproduces the real conditions of heat conduction in solids. Thus, about seventy years ago, Cattaneo proposed a modification of the Fourier equation now called the Cattaneo–Vernotte equation. This is the hyperbolic partial differential equation (PDE) and contains the parameter τq called the relaxation time (the lag time of the heat flux in relation to the temperature gradient) [1,2]. Especially important differences between the Fourier model and the real course of thermal processes appear in the case of microscale heat transfer problems. For example, the very high heating rates accompanying the heating of thin metal films with a laser beam mean that the inclusion of the finite value of thermal wave velocity must be taken into account. The deviations appear mainly when the mean free path of the heat carriers becomes comparable to the characteristic length of the domain considered and the time scale of interest becomes comparable to or smaller than the relaxation time of the heat carriers [3,4].
For the analysis of this type of process, the model with two delay times called a dual-phase lag model is presently applied. In addition to the relaxation time, the thermalization time is introduced. The relaxation time τq takes into account the small-scale response in time, while the thermalization time τT takes into account the small-scale response in space [3,4,5,6]. The dual phase lag equation (DPLE) results from the generalized form of the Fourier law
q ( x , t + τ q ) = λ   T ( x , t + τ T )
where q is a heat flux vector; ∇T is a temperature gradient; λ is a thermal conductivity; and x and t denote the geometrical co-ordinates and time, respectively. Both sides of the last dependence are developed into a power series and, finally (depending on the number of components), the first- or second-order DPLE can be obtained.
The literature on equations with two delay times is very extensive (especially in the case of the first order equations), and here we quote only a few important articles. The first publications concerning the model with two delay times appeared in the early nineties of the last century. There may be mentioned, for example, the papers [7,8,9]. Currently, one can already find books devoted to this type of non-Fourier heat conduction model, for example, [10,11,12,13].
In this brief literature review, the selected papers on analytical and numerical solutions of first-order DPLE will be listed. First of all, the works containing the analytical solutions of first-order DPLE (usually 1D tasks) will be mentioned [14,15,16,17,18]. The paper [14] concerns the laser heating of ultra-thin metal film; in the papers [15,16], the bio-heat transfer problems are discussed; while in the paper [18], the multi-layered metal domain is considered.
A much larger number of articles concern the application of numerical methods in the tasks based on the models with two delay times. It should be pointed out that, first of all, different variants of finite difference method (FDM) are applied (see, for example, [19,20,21,22,23,24]). In the paper [19], the numerical model of heating of the double-layered thin film has been applied for the analysis of the thermal deformation process. In the paper [20], the 3D FDM numerical model of the thin metal film heating has been presented. In [21], the explicit scheme of the FDM has been used and the problem of biological tissue freezing process has been discussed. The stability problem of the algorithm of this type is analyzed in [22]. In [23], the problem based on DPLE has been solved using the alternating direction implicit FDM scheme. In the paper [24], the adaptation of typical boundary conditions for non-Fourier equations has been shown. The FDM numerical solutions of the inverse problems can also be found (e.g., [25]).
The number of works presenting solutions using the other numerical methods is significantly smaller. Here, one can mention the papers [26,27,28,29,30,31,32]. In particular, solutions based on the finite element method [26,27,28], the boundary element method [29], the control volume method [30], or the lattice Boltzmann method [31,32] are discussed in the above-mentioned papers.
Literature on the second-order DPLE is not as extensive as for the first-order equations. As an example, the papers [33,34,35,36,37] can be mentioned. The main subject of these papers (except [37]) is related to the construction of algorithms for numerical modeling of problems described by second-order DPLE (the different variants of FDM are used). Additionally, the transformed second-order equations are shown in [20,34,37] and the changed forms are more convenient at the stage of numerical modeling.
At the stage of melting/resolidification modeling, the concept presented in [38] (for the macroscale problems) is applied. The capacity of the internal heat source related to the phase change is proportional to the melting/solidification rate. To define this function (in particular, the volumetric fraction of liquid state fL(T)) in the form of a continuous one, the melting point corresponding to the temperature Tm is conventionally replaced by a certain interval [Tm − ∆T, Tm + ∆T], and then the function discussed can be described by a broken line. For this interval, the substitute thermal capacity is defined and the one domain approach can be used. It should be pointed out that the testing computations show a little impact of the interval ∆T width (within reasonable limits) on the results of numerical simulations.
At the beginning of the part of the article devoted to own research, the assumed form of the dual- phase lag equation and the mathematical formulas determining the laser action and the evolution of phase change latent heat are presented. Both phenomena are taken into account by an introduction to the energy equation of the functions determining the efficiency of internal heat sources. Next, the numerical algorithm based on the implicit scheme of FDM is discussed. In the final part of the paper, the results of numerical computations concerning the heating/cooling process of the thin metal film made of chrome are shown. The conclusions resulting from the performed research are also formulated.

2. Governing Equations

The starting point for the formulation of the energy equation with delays is the generalized Fourier law (1). To obtain the DPLE, the left and right sides of Equation (1) are developed into the Taylor series
q ( x ,   t ) + τ q   q ( x ,   t )   t + τ q 2 2 2   q ( x ,   t )   t 2 + = λ [ T ( x ,   t ) + τ T     T ( x ,   t )   t + τ T 2 2 2 T ( x ,   t )   t 2 + ]
Let us apply the well-known diffusion equation, namely,
c   T ( x ,   t )   t =   ·   q   ( x ,   t ) + Q ( x , t )
where c is a volumetric specific heat and Q(x, t) is a capacity of volumetric internal heat sources.
When the components containing the second derivatives (Equation (2)) are taken into account, after mathematical manipulations, Equation (3) takes the form
x Ω :   c [   T ( x ,   t )   t + τ q 2   T ( x ,   t )   t 2 + τ q 2 2   3 T ( x ,   t )     t 3 ]   =   [ λ   T ( x ,   t ) ] + τ T   { [ λ   T ( x ,   t ) ] }   t + τ T 2 2   2 { [ λ   T ( x ,   t ) ] }     t 2   + Q ( x , t ) + τ q     Q ( x ,   t )   t + τ q 2 2 2   Q ( x ,   t )   t 2  
When the melting and resolidification problem is considered, the internal heat source in Equation (4) must contain the term controlling the phase change process. This appropriate source function Qm (x, t) can be defined as
Q m ( x , t ) = L   f L ( x ,   t )   t L τ q 2 f L ( x ,   t )   t 2 L τ q 2 2 3 f L ( x , t )   t 3
where L is the volumetric latent heat phase change and fL(x, t) is the volumetric molten state fraction at the neighborhood of the point considered. The last equation is a generalization of what is well known in the thermal theory of foundry processes definition of Qm (e.g., [38]). The function fL(x, t) is equal to zero at the beginning of the heating process until T1 = Tm − ∆T and fL(x, t) = 1 for T2 > Tm + ∆T (Tm is the melting point). In the interval [Tm − ∆T, Tm + ∆T], the function fL(x, t) changes from 0 to 1 in a linear way (such an assumption is fully acceptable). Generally speaking, the volumetric liquid state fraction is given in the form of broken line, this means
f L ( x , t ) = { 1 T ( x , t ) > T m + Δ T T ( x , t ) T m + Δ T 2 Δ T T m Δ T T ( x , t ) T m + Δ T 0 T ( x , t ) < T m Δ T
The derivative of fL(x, t) with respect to the temperature is equal to 0 for T(x, t) < Tm − ∆T and T(x, t) > Tm +T, while between the border temperatures dfL(x, t)/dT = 1/2∆T. Thus, the source term Qm(x, t) acts only for T(x, t) from the interval [Tm − ∆T, Tm + ∆T], and then
Q m ( x , t ) = L 2 Δ   T [   T ( x ,   t )   t + τ q   T 2 ( x ,   t )   t 2 + τ q 2 2   T 3 ( x ,   t )   t 3 ] ,     T ( x ,   t ) [ T m Δ T , T m + Δ T ]
Let us introduce the piece-vise constant function C(T)
C ( T ) = { c 2 T ( x , t ) > T m + Δ T 0.5 ( c 1 + c 2 ) + L 2 Δ T T m Δ T T ( x , t ) T m + Δ T c 1 T ( x , t ) < T m Δ T
where c1 and c2 are the volumetric specific heats of the solid and liquid states, respectively. Then, Equation (4) can be written as follows
x Ω :   C ( T ) [   T ( x ,   t )   t + τ q 2   T ( x ,   t )   t 2 + τ q 2 2   3 T ( x ,   t )     t 3 ]   =   [ λ   T ( x ,   t ) ] + τ T   { [ λ   T ( x ,   t ) ] }   t + τ T 2 2   2 { [ λ   T ( x ,   t ) ] }     t 2   + Q l ( x , t ) + τ q     Q l ( x ,   t )   t + τ q 2 2 2   Q l ( x ,   t )   t 2  
Thermal conductivity λ in Equation (9) is defined just like the parameter C(T).
The mathematical formula determining the intensity of the internal heat source Ql (x, t) resulting from the laser action can be taken in the form [39]
Q l ( x ,   t ) = ( 1 R )   I 0 δ   t p   exp [ x 1 2 + x 2 2 r D 2   x 3 δ 4 ln 2 ( t 2 t p ) 2 t p 2 ] ,   x = { x 1 , x 2 , x 3 }
where I0 [J/m2] is the laser intensity, tp [s] is the characteristic time of laser pulse, δ [m] is the optical penetration depth, R is the reflectivity of the irradiated surface, rD [m] is the laser beam radius, and x3 is a vertical axis. The derivatives of Ql with respect to time can be found analytically.
On the outer surface of the system, the adiabatic conditions are assumed (the external heat flux is taken into account in the appropriate source function). The mathematical form of the Neumann boundary condition for the second-order DPLE is as follows [36]
x Γ :    λ [ n · T ( x ,   t ) + τ T [ n · T ( x ,   t ) ] t + τ T 2 2 2 [ n · T ( x ,   t ) ] t 2 ] = q b ( x ,   t ) + τ q q b ( x ,   t ) t + τ q 2 2 2   q b ( x ,   t )   t 2
where n is a normal outward vector and (in the case considered) qb (x, t) = 0, of course.
The mathematical model is also supplemented by the initial conditions
t = 0 :     T ( x , 0 ) = T p ,     T ( x , t ) t | t = 0 = Q ( x , 0 ) c 1 ,   2 T ( x , t ) t 2 | t = 0 = 1 c 1 Q ( x , t ) t | t = 0
where Tp is an initial temperature.

3. Mathematical Description of 1D Problem

At the stage of numerical modeling, the 1D problem was considered and the basis for the construction of the FDM algorithm is the following system of equations and conditions:
-
energy equation for thin metal film domain
C ( T ) [   T ( x ,   t )   t + τ q 2   T x , t )   t 2 + τ q 2 2   3 T ( x ,   t )     t 3 ] =     x [ λ   T ( x ,   t )   x ] +   τ T 2     t     x [ λ   T ( x ,   t )   x ] + τ T 2 2   3     t 2     x [ λ   T ( x ,   t )   x ]   + Z ( x , t )
where
Z ( x , t ) = Q l ( x , t ) + τ q   Q l    ( x ,   t )   t + τ q 2 2 2   Q l ( x , t )   t 2
-
source function Ql
Q l ( x ,   t ) = ( 1 R )   I 0 δ   t p   exp [   x δ 4 ln 2 ( t 2 t p ) 2 t p 2 ]
adiabatic boundary conditions
x = 0 G :        T ( x ,   t ) + τ T [   T ( x ,   t )   x ] t + τ T 2 2 2 [   T ( x ,   t )   x ] t 2 = 0
-
initial conditions
t = 0 :     T ( x , 0 ) = T p ,     T ( x , t ) t | t = 0 = Q ( x , 0 ) c 1 ,   2 T ( x , t ) t 2 | t = 0 = 1 c 1 Q ( x , t ) t | t = 0

4. Numerical Model Based on FDM

Let T i f = T ( x i , f Δ t ) , where x i = i h , i = 0, 1, …, n and f = 0, 1,…, F. Taking into account the initial conditions (17), one has
T i 0 = T p T i 1 T i 0 Δ t = 1 c 1 Q l ( x i , 0 )     T i 1 = T i 0 + Δ t c 1 Q l ( x i , 0 )   T i 2 2 T i 1 + T i 0 ( Δ t ) 2 = 1 c 1 Q l ( x i , t ) t | t = 0     T i 2 = 2 T i 1 T i 0 + ( Δ t ) 2 c 1 Q j ( x i , t ) t | t = 0
For the transition t f−1t f (f ≥ 3), the approximate form of Equation (13) resulting from the introduction of the assumed differential quotients is as follows
C i f 1 [   T i f T i f 1 Δ   t + τ q T i f 2 T i f 1 + T i f 2 ( Δ t ) 2 + τ q 2 2 T i f 3 T i f 1 + 3 T i f 2 T i f 3 ( Δ t ) 3 ] = [     x ( λ   T   x ) ]   i   f + τ T Δ   t { [     x ( λ   T   x ) ]   i   f [     x ( λ   T   x ) ]   i   f 1 } + w T τ T 2 2 ( Δ t ) 2   { [     x ( λ   T   x ) ]   i   f 2 [     x ( λ   T   x ) ]   i   f 1 + [     x ( λ   T   x ) ]   i   f 2 }   + Z i f
where
[     x ( λ   T   x ) ] i s = 1 h [ ( λ   T   x ) i + 0.5 s ( λ   T   x ) i 0.5 s ] = 1 h ( λ i + 0.5 s T i + 1   s T i   s h λ i 0.5 s T i   s T i 1   s h ) = 1 h ( λ i s + λ i + 1 s 2 T i + 1   s T i   s h λ i 1 s + λ i s 2 T i   s T i 1   s h ) = 1 2 h 2 [ ( λ i s + λ i + 1 s ) ( T i + 1   s T i   s ) + ( λ i 1 s + λ i s ) ( T i 1   s T i   s ) ]
while s = f, s = f − 1, or s = f − 2.
Thus,
  T i f T i f 1 Δ   t + τ q T i f 2 T i f 1 + T i f 2 ( Δ t ) 2 + τ q 2 2 T i f 3 T i f 1 + 3 T i f 2 T i f 3 ( Δ t ) 3 = 2 ( Δ t ) 2 + 2 τ T Δ   t + τ T 2 4 C i f 1 h 2 ( Δ t ) 2 [ ( λ i f 1 + λ i + 1 f 1 ) ( T i + 1   f T i   f ) + ( λ i 1 f 1 + λ i f 1 ) ( T i 1   f T i   f ) ] 2 τ T Δ   t + 2 τ T 2 4 C i f 1 h 2 ( Δ t ) 2 [ ( λ i f 1 + λ i + 1 f 1 ) ( T i + 1   f 1 T i   f 1 ) + ( λ i 1 f 1 + λ i f 1 ) ( T i 1   f 1 T i   f 1 ) ] + τ T 2 4 C i f 1 h 2 ( Δ t ) 2   [ ( λ i f 2 + λ i + 1 f 2 ) ( T i + 1   f 2 T i   f 2 ) + ( λ i 1 f 2 + λ i f 2 ) ( T i 1   f 2 T i   f 2 ) ] + Z i f C i f 1
and after mathematical manipulations, one has
A i f T i 1 f + B i f T i f + C i f T i + 1 f = D i f
where
A i f = [ 2 ( Δ t ) 2 + 2 τ T Δ   t + τ T 2 ] ( λ i 1 f 1 + λ i f 1 ) 4 C i f 1 h 2 ( Δ t ) 2 , B i f = 2 ( Δ t ) 2 + 2 τ q Δ   t + τ q 2 2 ( Δ t ) 3 + [ 2 ( Δ t ) 2 + 2 τ T Δ   t + τ T 2 ] ( λ i 1 f 1 + 2 λ i f 1 + λ i + 1 f 1 ) 4 C i f 1 h 2 ( Δ t ) 2 , C i f = [ 2 ( Δ t ) 2 + 2 τ T Δ   t + τ T 2 ] ( λ i + 1 f 1 + λ i f 1 ) 4 C i f 1 h 2 ( Δ t ) 2 , D i f = 2 ( Δ t ) 2 + 4 τ q   Δ t + 3 τ q 2 2 ( Δ t ) 3 T i f 1 2 τ q   Δ t + 3 τ q 2 2 ( Δ t ) 3 T i f 2 + τ q 2 2 ( Δ t ) 3 T i f 3 2 τ T Δ t + 2 τ T 2 4 C i f 1 h 2 ( Δ t ) 2 [ ( λ i f 1 + λ i + 1 f 1 ) ( T i + 1   f 1 T i   f 1 ) + ( λ i 1 f 1 + λ i f 1 ) ( T i 1   f 1 T i   f 1 ) ] + τ T 2 4 C i f 1 h 2 ( Δ t ) 2   [ ( λ i f 2 + λ i + 1 f 2 ) ( T i + 1   f 2 T i   f 2 ) + ( λ i 1 f 2 + λ i f 2 ) ( T i 1   f 2 T i   f 2 ) ] + Z i f C i f 1
The approximation of boundary conditions (16) is as follows
-
for x = 0:
T 1 f T 0 f h + τ T Δ t ( T 1 f T 0 f h T 1 f 1 T 0 f 1 h ) + τ T 2 2 ( Δ t ) 2 ( T 1 f T 0 f h 2 T 1 f 1 T 0 f 1 h + T 1 f 2 T 0 f 2 h ) = 0
-
for x = G:
T n f T n 1 f h + τ T Δ t ( T n f T n 1 f h T n f 1 T n 1 f 1 h ) + τ T 2 2 ( Δ t ) 2 ( T n f T n 1 f h 2 T n f 1 T n 1 f 1 h + T n f 2 T n 1 f 2 h ) = 0
or
[ 2 ( Δ t ) 2 + 2 τ T Δ t + τ T 2 ] T 0 f + [ 2 ( Δ t ) 2 + 2 τ T Δ t + τ T 2 ] T 1 f = ( 2 τ T Δ t + 2 τ T 2 ) ( T 1 f 1 T 0 f 1 ) τ T 2 ( T 1 f 2 T 0 f 2 )
and
[ 2 ( Δ t ) 2 + 2 τ T Δ t + τ T 2 ] T n 1 f + [ 2 ( Δ t ) 2 + 2 τ T Δ t + τ T 2 ] T n f = ( 2 τ T Δ t + 2 τ T 2 ) ( T n f 1 T n 1 f 1 ) τ T 2 ( T n f 2 T n 1 f 2 )
The transition from t f − 1 to t f (f ≥ 3) requires the solution of the system of Equations (22), (26), and (27) with a three-band main matrix that can be solved quickest using the Thomas algorithm [40].

5. Results of Computations

The computations were performed for the thin metal film of the thickness G = 100 nm made of chromium. The surface x = 0 is subjected to the laser pulse with the parameters R = 0.93, I0 = 3825 J/m2, δ = 15.3 nm, and tp = 10 ps—c.f. Formula (15). The following values of chromium parameters are assumed: thermal conductivities λ1 = 93 W/(m K), λ2 = 35 W/(m K), volumetric specific heats c1 = 3.2148 MJ/(m3 K), c2 = 2.79276 MJ/(m3 K) [24], relaxation time τq = 0.136 ps, thermalization time τT = 7.86 ps, melting temperature Tm = 2180 K, and volumetric heat of fusion L = 2904 MJ/m3.
Different temporal-spatial meshes were considered; see Table 1. For each combination of mesh steps, the temperature after 80 ps on the heated surface was recorded. As can be seen, the result does not depend much on the time step, while the number of nodes is important. Analysis of the results obtained showed that the satisfactory accuracy provides the geometrical mesh containing n = 1000 nodes (h = 0.1 nm) and time step ∆t = 0.0005 ps (see Table 1, column 5).
Figure 1 shows the temperature courses at the points x = 0 (heated surface) and x = 20 nm. The results were compared to those published in [24], where the first order DPL equation and slightly different partial melting model were used. As one can see, the differences are small, but visible. For the second-order model (solid lines), the maximum temperature of the heated surface occurs after 28.8 ps and is equal to 2323.28 K, while for the first-order model (dashed lines), the maximum temperature occurs after 28 ps and equals 2307.20 K.
In Figure 2, the fragments of heating/cooling curves for different widths of ΔT are shown. One can see that the results are similar. For ΔT = 3 K, the maximum temperature is the highest; for ΔT = 7 K, the maximum temperature is the lowest; while the biggest differences between them are of the order 10 K. After the resolidification process, the cooling curves for all cases are very similar. Further computations were carried out for ΔT = 3 K.
In Figure 3, the temperature courses on the irradiated surface for the laser intensity I0 = 3825 J/m2 and different characteristic times of laser pulse tp are presented. The impact of changes in this parameter on the course of the heating/cooling curves is clearly visible. It is obvious that an increase in intensity of the laser pulse causes a more rapid heating process. Interesting, however, are the effects of changes of characteristic time tp. Shortening this time increases the heating rate of the metal film and maxima of the cooling/heating curves shift to the left. At the same time, the maximum values of temperatures increase slightly. A detailed analysis of Equation (10) determining the time-dependent efficiency of the laser-related internal heat source shows that such an effect could be expected.
Finally, Figure 4 shows the changes of molten layer thicknesses in time. The successive curves correspond to tp = 2, 6, and 10 ps. For tp = 2 ps, the maximal depth is equal to g = 17.6 nm; for tp = 6 ps, it is equal to g = 15.3 nm; while for tp = 10 ps, it is equal to g = 13 nm.
Changes in the characteristic time of laser pulse cause a similar effect to that seen in the previous figure. The depth of the molten layer increases with the reduction of time tp and the process is more dynamic. Such analysis and information obtained on its basis can often be useful in engineering practice.

6. Conclusions

The mathematical model, the numerical algorithm, and exemplary results of the computations concerning the heating of thin metal film subjected to the laser beam are discussed, wherein the problem is described using the second-order DPLE. The laser beam power is so high that, in the domain under consideration, the partial melting and then resolidification of the material take place. To our knowledge, the application of the second-order DPLE to solve this type of the problem has not yet been presented. The task is treated as a non-linear one and the thermophysical parameters of the material are temperature-dependent. The melting point is substituted by a certain interval of temperature. This assumption allows one to introduce the continuous function determining the local and temporary volumetric liquid state fraction of the metal. The influence of the width of this interval on the results of numerical calculations is also examined. The laser action and the evolution of latent heat are taken into account by the introduction of two internal heat sources to the energy equation.
At the stage of numerical modeling, the 1D problem was considered, but the generalization of the algorithm and computer program implementing the numerical computations for the larger number of dimensions is not difficult. The computer program is based on the implicit scheme of the finite difference method. The algorithm is unconditionally stable and the transition from time t to t + ∆t requires the solution of the system of linear equations, whose main matrix is a three-diagonal one. In this place, the very effective Thomas algorithm was used.
The obtained results are in line with the expectations; their comparison with the solution of a similar problem described by the first-order equation shows visible, but slight differences.
We are planning the following further research in this scope:
-
numerical algorithm and computer program for 3D problems;
-
numerical algorithm and computer program for axially-symmetrical tasks (this geometry is very convenient because of the typical shape of the function describing the laser action);
-
adaptation of the algorithm presented in this work for modeling the ablation process;
-
research on other approaches to phase changes modeling.

Author Contributions

Conceptualization, E.M. and B.M.; formal analysis, E.M. and B.M.; software, E.M.; investigation, E.M. and B.M., writing—original draft, B.M. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Cattaneo, M.C. Sulla Conduzione de Calore. Atti de Seminario Matematico e Fisico; Della Universita di Modena: Modena, Italy, 1948; pp. 3–21. [Google Scholar]
  2. Cattaneo, M.C. A form of heat conduction equation which eliminates the paradox of instantaneous propagation. Comptes Rendus 1958, 247, 431–433. [Google Scholar]
  3. Escobar, R.A.; Ghai, S.S.; Jhon, M.S.; Amon, C.H. Multi-length and time scale thermal transport using the l lattice Boltzmann method with applications to electronics cooling. Int. J. Heat Mass Transf. 2006, 49, 97–107. [Google Scholar] [CrossRef]
  4. Flik, M.I.; Choi, B.I.; Goodson, K.E. Heat transfer regimes in microstructures. J. Heat Transf. 1992, 114, 664–674. [Google Scholar] [CrossRef] [Green Version]
  5. Tzou, T. A unified field approach for heat conduction from macro- to micro-scales. J. Heat Transf. 1995, 117, 8–16. [Google Scholar] [CrossRef]
  6. Cabrera, J.; Castro, M.A.; Rodríquez, F.; Martín, J.A. Difference schemes for numerical solutions of lagging models of heat conductions. Math. Comupt. Model. 2013, 57, 1625–1632. [Google Scholar] [CrossRef]
  7. Özişik, M.N.; Tzou, D.Y. On the wave theory in heat conduction. J. Heat Transf. 1994, 116, 526–535. [Google Scholar] [CrossRef]
  8. Orlande, R.B.; Özişik, M.N.; Tzou, D.Y. Inverse analysis for estimating the electron-phonon coupling factor in thin metal films. J. Appl. Phys. 1995, 78, 1843–1849. [Google Scholar] [CrossRef]
  9. Al Nimr, M.A. Heat transfer mechanisms during short duration laser heating of thin metal films. Int. J. Thermophys. 1997, 18, 1257–1268. [Google Scholar] [CrossRef]
  10. Zhang, Z.M. Nano/Microscale Heat Transfer; McGraw-Hill: New York, NY, USA, 2007. [Google Scholar]
  11. Tzou, D.Y. Macro- to Microscale Heat Transfer. The Lagging Behavior; John Wiley & Sons: Hoboken, NJ, USA, 2015. [Google Scholar]
  12. Faghri, A.; Zhang, Y.; Howell, J. Advanced Heat and Mass Transfer; Global Digital Press: Columbia, MO, USA, 2010. [Google Scholar]
  13. Smith, A.N.; Norris, P.M. Microscale Heat Transfer, Chapter 18. In Heat Transfer Handbook; John Willey & Sons: Hoboken, NJ, USA, 2003. [Google Scholar]
  14. Ciesielski, M. Analytical solution of the dual phase lag equation describing the laser heating of thin metal film. J. Appl. Math. Comput. Mech. 2017, 16, 33–40. [Google Scholar] [CrossRef] [Green Version]
  15. Kumar, S.; Srivastava, A. Finite integral transform-based analytical solutions of dual phase lag bio-heat transfer equation. Appl. Math. Model. 2017, 52, 378–403. [Google Scholar] [CrossRef]
  16. Liu, K.-C.; Wang, J.-C. Analysis of thermal damage to laser irradiated tissue based on the dual-phase-lag model. Int. J. Heat Mass Transf. 2014, 70, 621–628. [Google Scholar] [CrossRef]
  17. Mohammadi-Fakhar, V.; Momeni-Masuleh, S.H. An approximate analytic solution of the heat conduction equation at nanoscale. Phys. Lett. A 2010, 374, 595–604. [Google Scholar] [CrossRef]
  18. Ramadan, K. Semi-analytical solutions for the dual phase lag heat conduction in multi-layered media. Int. J. Therm. Sci. 2009, 48, 14–25. [Google Scholar] [CrossRef]
  19. Wang, H.; Dai, W.; Melnik, R. A finite difference method for studying thermal deformation in a double-layered thin film exposed to ultrashort pulsed lasers. Int. J. Therm. Sci. 2006, 45, 1179–1196. [Google Scholar] [CrossRef]
  20. Dai, W.; Nassar, R. A compact finite difference scheme for solving a three-dimensional heat transport equation in a thin film. Numer. Methods Partial Differ. Equ. 2000, 16, 441–458. [Google Scholar] [CrossRef]
  21. Mochnacki, B.; Majchrzak, E. Numerical model of thermal interactions between cylindrical cryoprobe and biological tissue using the dual-phase lag equation. Int. J. Heat Mass Transf. 2017, 108, 1–10. [Google Scholar] [CrossRef]
  22. Majchrzak, E.; Mochnacki, B. Dual-phase-lag equation. Stability conditions of a numerical algorithm based on the explicit scheme of the finite difference method. J. Appl. Math. Comput. Mech. 2016, 15, 89–96. [Google Scholar] [CrossRef] [Green Version]
  23. Ciesielski, M. Application of the alternating direction implicit method for numerical solution of the dual phase lag equation. J. Theor. Appl. Mech. Pol. 2017, 55, 839–852. [Google Scholar] [CrossRef] [Green Version]
  24. Majchrzak, E.; Mochnacki, B. Dual-phase lag model of thermal processes in a multi-layered microdomain subjected to a strong laser pulse using the implicit scheme of FDM. Int. J. Therm. Sci. 2018, 133, 240–251. [Google Scholar] [CrossRef]
  25. Mochnacki, B.; Paruch, M. Estimation of relaxation and thermalization times in microscale heat transfer model. J. Theor. Appl. Mech. Pol. 2013, 51, 837–845. [Google Scholar]
  26. Kumar, P.; Kumar, D.; Rai, K.N. A numerical study of dual-phase-lag model of bio-heat transfer during hyperthermia treatment. J. Therm. Biol. 2015, 49–50, 98–105. [Google Scholar] [CrossRef] [PubMed]
  27. Kumar, D.; Kumar, P.; Rai, K.N. A study on DPL model of heat transfer in bi-layer tissues during MFH Treatment. Comput. Biol. Med. 2016, 75, 160–172. [Google Scholar] [CrossRef] [PubMed]
  28. Kumar, D.; Rai, K.N. A study on thermal damage during hyperthermia treatment based on DPL model for multilayer tissues using finite element Legendre wavelet Galerkin approach. J. Therm. Biol. 2016, 62, 170–180. [Google Scholar] [CrossRef] [PubMed]
  29. Majchrzak, E. Numerical solution of dual phase lag model of bioheat transfer using the general boundary element method. CMES Comput. Model. Eng. Sci. 2010, 69, 43–60. [Google Scholar]
  30. Mochnacki, B.; Ciesielski, M. Micro-scale heat transfer. Algorithm basing on the Control Volume Method and the identification of relaxation and thermalization times using the search method. CMMS 2015, 15, 353–361. [Google Scholar]
  31. Patidar, S.; Kumar, S.; Srivastava, A.; Singh, S. Dual phase lag model-based thermal analysis of tissue phantoms using lattice Boltzmann method. Int. J. Therm. Sci. 2016, 103, 41–56. [Google Scholar] [CrossRef]
  32. Ho, J.R.; Kuo, C.P.; Jiaung, W.S. Study of heat transfer in multilayered structure within the framework of dual-phase-lag heat conduction model using lattice Boltzmann method. Int. J. Heat Mass Transf. 2003, 46, 55–69. [Google Scholar] [CrossRef]
  33. Castro, M.A.; Rodríques, F.; Cabrera, J.; Martín, J.A. A compact difference scheme for numerical solution of second order dual-phase-lagging models of microscale heat transfer. J. Comput. Appl. Math. 2016, 291, 432–440. [Google Scholar] [CrossRef] [Green Version]
  34. Deng, D.; Jiang, Y.; Liang, D.L. High-order finite difference methods fora second order dual-phase-lagging models of microscale heat transfer. Appl. Math. Comput. 2017, 309, 31–48. [Google Scholar]
  35. Majchrzak, E.; Mochnacki, B. Implicit scheme of the finite difference method for a second-order dual phase lag equation. J. Theor. Appl. Mech. Pol. 2018, 56, 393–402. [Google Scholar] [CrossRef] [Green Version]
  36. Majchrzak, E.; Mochnacki, B. Numerical solutions of the second-order dual-phase lag equation using the explicit and implicit schemes of the finite difference method. Int. J. Numer. Methods Heat Fluid Flow 2019, 20, 2099–2120. [Google Scholar] [CrossRef]
  37. Askarizadeh, H.; Baniasadi, E.; Ahmadikia, H. Equilibrium and nonequilibrium thermodynamic analysis of high-order dual-phase-lag heat conduction. Int. J. Heat Mass Transf. 2017, 104, 301–309. [Google Scholar] [CrossRef]
  38. Ciesielski, M.; Mochnacki, B. Comparison of approaches to the numerical modelling of pure metals solidification using the control volume method. Int. J. Cast Metals Res. 2019, 32, 213–220. [Google Scholar] [CrossRef]
  39. Grigoropoulos, C.P.; Chimmalgi, A.; Hwang, D.J. Laser Ablation and Its Applications; Springer Series in Optical Sciences; Springer: New York, NY, USA, 2007; pp. 473–504. [Google Scholar]
  40. Datta, B.N. Numerical Linear Algebra and Applications, 2nd ed.; SIAM: Philadelphia, PA, USA, 2010; p. 162. [Google Scholar]
Figure 1. Temperature histories at the points x = 0 and x = 20 nm; comparison of the results obtained using the first-order (dashed line [24]) and second-order dual phase lag (DPL) (solid line) for I0 = 3825 J/m2, tp = 10 ps.
Figure 1. Temperature histories at the points x = 0 and x = 20 nm; comparison of the results obtained using the first-order (dashed line [24]) and second-order dual phase lag (DPL) (solid line) for I0 = 3825 J/m2, tp = 10 ps.
Mathematics 08 00999 g001
Figure 2. Comparison of obtained temperature courses for different intervals ΔT (on the surface x = 0), (1) ΔT = 3 K, (2) ΔT = 5 K, (3) ΔT = 7 K.
Figure 2. Comparison of obtained temperature courses for different intervals ΔT (on the surface x = 0), (1) ΔT = 3 K, (2) ΔT = 5 K, (3) ΔT = 7 K.
Mathematics 08 00999 g002
Figure 3. Heating/cooling curves on the surface x = 0 for different characteristic times of laser pulse. tp: (1) tp = 2 ps, (2) tp = 6 ps, (3) tp = 10 ps (I0 = 4500 J/m2).
Figure 3. Heating/cooling curves on the surface x = 0 for different characteristic times of laser pulse. tp: (1) tp = 2 ps, (2) tp = 6 ps, (3) tp = 10 ps (I0 = 4500 J/m2).
Mathematics 08 00999 g003
Figure 4. Changes in the thickness of the molten layer for different characteristic times of laser pulse, (1) tp = 2 ps, (2) tp = 6 ps, (3) tp = 10 ps.
Figure 4. Changes in the thickness of the molten layer for different characteristic times of laser pulse, (1) tp = 2 ps, (2) tp = 6 ps, (3) tp = 10 ps.
Mathematics 08 00999 g004
Table 1. Temperature of the irradiated surface after 80 ps for different time steps and number of nodes.
Table 1. Temperature of the irradiated surface after 80 ps for different time steps and number of nodes.
t [ps]n = 100n = 200n = 500n = 1000
0.00011484.221498.471507.061511.45
0.000251484.191498.431507.071511.43
0.00051484.081498.441507.051511.47

Share and Cite

MDPI and ACS Style

Majchrzak, E.; Mochnacki, B. Second-Order Dual Phase Lag Equation. Modeling of Melting and Resolidification of Thin Metal Film Subjected to a Laser Pulse. Mathematics 2020, 8, 999. https://doi.org/10.3390/math8060999

AMA Style

Majchrzak E, Mochnacki B. Second-Order Dual Phase Lag Equation. Modeling of Melting and Resolidification of Thin Metal Film Subjected to a Laser Pulse. Mathematics. 2020; 8(6):999. https://doi.org/10.3390/math8060999

Chicago/Turabian Style

Majchrzak, Ewa, and Bohdan Mochnacki. 2020. "Second-Order Dual Phase Lag Equation. Modeling of Melting and Resolidification of Thin Metal Film Subjected to a Laser Pulse" Mathematics 8, no. 6: 999. https://doi.org/10.3390/math8060999

APA Style

Majchrzak, E., & Mochnacki, B. (2020). Second-Order Dual Phase Lag Equation. Modeling of Melting and Resolidification of Thin Metal Film Subjected to a Laser Pulse. Mathematics, 8(6), 999. https://doi.org/10.3390/math8060999

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop