Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 55, 3, pp. 839-852, Warsaw 2017 DOI: 10.15632/jtam-pl.55.3.839 APPLICATION OF THE ALTERNATING DIRECTION IMPLICIT METHOD FOR NUMERICAL SOLUTION OF THE DUAL PHASE LAG EQUATION Mariusz Ciesielski Czestochowa University of Technology, Institute of Computer and Information Sciences, Częstochowa, Poland e-mail: mariusz.ciesielski@icis.pcz.pl The problem discussed in the paper concerns the numerical modeling of thermal processes proceeding in micro-scale described using the Dual Phase Lag Model (DPLM) in which the relaxation and thermalization time appear. The cylindrical domain of a thin metal film subjected to a strong laserpulsebeam is considered.The laser action is taken intoaccountby the introduction of an internal heat source in the energy equation.At the stage of numerical modeling, the ControlVolumeMethod is used and adapted to resolve the hyperbolic partial differential equation. In particular, the Alternating Direction Implicit (ADI) method for DPLM is presented and discussed. The examples of computations are also presented. Keywords:micro-scale heat transfer, dual phase lagmodel, control volumemethod, alterna- ting direction implicit method, numerical simulation 1. Introduction Thermal processes proceeding in the microscale are characterized, as a rule, by an extremely short duration, extreme temperature gradients and very small geometrical dimensions of the domain considered. It is a reason that typical mathematical models basing on the macroscopic Fourier-type equations are not suitable for the analysis of this type problems. In the recent years, the problem of heat transfer through domains subjected to an strong external heat source (e.g. an ultrafast laser pulse) has been of vital importance inmicrotechnology applications, and it is a reason that the problems connected with fast heating of solids has become a very active research area (Tzou, 2015; Zhang, 2007; Chen et al., 2004). From the mathematical point of view, nowadays there exist different models describing the mechanism of the process discussed. In this group, the microscopic two-temperature parabolic or hyperbolic models (belonging to a group of continuum models) should be mentioned (Chen andBeraun, 2001; Kaba andDai, 2005; Lin and Zhigilei, 2008;Majchrzak, 2012;Majchrzak and Dziatkiewicz, 2015). The two-step parabolic andhyperbolicmodels involve two energy equations determining the thermal processes in the electron gas and themetal lattice. The coupling factor combining these equations is introduced. Depending on the variant of the model, parabolic or hyperbolic PDEs are considered. Assuming certain simplifications, the two-temperature model can be transformed into a single equation containing the second order time derivative andhigher ordermixed derivative in both time and space (known as the dual phase lagmodel (DPLM)). In this equation, two positive constants τq, τT appear and they correspond to the relaxation time, which is the mean time for electrons to change their energy states and the thermalization time, which is the mean time required for electrons and lattice to reach equilibrium (Orlande et al., 1995). TheCattaneo-Vernotte and the dual phase lagmodels belong also to the group of continuum ones. They result from the generalization of the well-known Fourier law. To take into account the finite velocity of a thermalwave the lag timebetween the heat fluxand temperature gradient 840 M. Ciesielski has been introduced (Cattaneo, 1958). The Cattaneo-Vernotte hyperbolic equation (CVE) can be treated as a certain microscale heat transfer model, but for this purpose is rarely used. The model discussed often finds application in the case of bioheat transfer problems, e,g. (Ciesielski et al., 2016). In fact, according to literature, e.g. (Mitra et al., 1995) the lag time (relaxation time) for processedmeat is of the order of several seconds. Introduction of two lag times in the generalized form of the Fourier law (relaxation and thermalization ones) leads, after relatively simplemathematicalmanipulations, to thedual phase lag equation. At present, in literature one can find big number of analytical and (first of all) numerical solutions of various thermal problems described by this model. The majority of the solutions presented in the literature concerns the 1D problems. Such an assumption is often fully acceptable. For example, considering the laser pulse interactions with thin metal films it is reasonable to treat the interactions as a one-dimensional heat transfer process (Chen and Beraun, 2001). In this paper, the axially-symmetrical problem is analyzed. Most of theworks in thisfield concernsdirectproblems.Homogeneous andalsoheterogeneous domains are considered. The problem of the single layer heating was discussed, among others, by Tang and Araki (1999), Kaba and Dai (2005), Mochnacki and Ciesielski (2012), Majchrzak and Turchan (2016). In the subject of non-homogeneous micro-domains, one can mention the paper presented by Dai and Nassar (2000), in which the heat transfer in a double layered gold- chromium film is analyzed, and the papers prepared by Majchrzak et al. (2009a,b) concerning a multi-layered film subjected to ultrafast laser heating. Both in the case of the CVE andDPLE, the typical boundary conditions appearing in heat transfer problems should bemodified in a adequate way. In literature, one can find works devoted to sensitivity of the transient temperature field in microdomainswith respect to the dual phase lagmodel parameters (Majchrzak, andMochnacki, 2014). The issue of the inverse problems was also developed, e.g. by Mochnacki and Paruch (2013), Dziatkiewicz et al. (2014), Mochnacki and Ciesielski (2015). A next group of microscale heat transfer models is based on the Boltzmann transport equ- ation (BTE). It is a conservation equationwhere theconservedquantity is thenumberofparticles (Tian and Yang, 2008). The general form of BTE is rather complex, but it can be modified to analyze special tasks, for instance systems created by phonons, electrons, etc. In this field, de- serving special attention is repeatedly cited paper presented by Escobar et al. (2006). One can also mention the work by Belhayat-Piasecka and Korczak (2016) in which the microscale heat transport was analyzed using the interval lattice Boltzmann method. Microscale heat transfer processes can be also considered using the molecular approaches (Smith and Norris, 2003; Theodosiou and Saravanos, 2007; Chen et al., 2007; Liu and Tsai, 2009). 2. Governing equations Let us consider the diffusion equation in the domainΩ (r,z)∈Ω c∂T(r,z,t) ∂t =−∇·q(r,z,t)+Q(r,z,t) (2.1) where c= c(T) is the volumetric specific heat, q(r,z,t) is the heat flux vector, Q(r,z,t) is the capacity of internal heat sources, r, z, t are the geometrical co-ordinates and time. The value of heat flux is determined by Tzou’s dual-phase-lag theory (Tzou, 2015), as the generalization of the Fourier law, in particular q(r,z,t+ τq)=−λ∇T(r,z,t+ τT) (2.2) Application of the alternating direction implicit method... 841 where τq is called the relaxation time, while τT is the thermalization time, λ = λ(T) is the thermal conductivity, ∇T(r,z,t) is the temperature gradient. For τT = 0, this model leads to the Cattaneo-Vernotte equation, while for τT =0 and τq =0 it corresponds to the Fourier law. The Taylor series expansions of equation (2.2) is the following q(r,z,t)+ τq ∂q(r,z,t) ∂t =−λ [ ∇T(r,z,t)+ τT ∂∇T(r,z,t) ∂t ] (2.3) Introducing formula (2.3) into equation (2.1) one obtains c [∂T(r,z,t) ∂t + τq ∂2T(r,z,t) ∂t2 ] =∇· [λ∇T(r,z,t)]+ τT ∂∇· [λ∇T(r,z,t)] ∂t +Q(r,z,t)+ τq ∂Q(r,z,t) ∂t (2.4) In the case of theaxially-symmetrical taskdiscussed in thiswork, the component∇·[λ∇T(r,z,t)] is the following ∇· [λ∇T(r,z,t)] = 1 r ∂ ∂r [ rλ ∂T(r,z,t) ∂r ] + ∂ ∂z [ λ ∂T(r,z,t) ∂z ] (2.5) It should be pointed out that the boundary conditions (which appear in the typical Fourier heat conduction models) for the DPL should be transformed to the form (r,z)∈Γ : qb(r,z,t)+ τq ∂qb(r,z,t) ∂t =−λ [ n ·∇T(r,z,t)+ τT ∂[n ·∇T(r,z,t)] ∂t ] (2.6) In Fig. 1, the domain considered (limited by the planes z = 0, z = Z and surface r = R) is shown. Fig. 1. Cylindrical micro-domain The effects of femtosecond laser pulse irradiation on the upper surface limiting the system causes that the energy is delivered tometal and its absorption occurs. The internal heat source Q(r,z,t) generated inside metal is related with action of the laser beam (Chen and Beraun, 2001) Q(r,z,t) = √ β π 1−Rf tpδ I0exp ( −z δ ) exp ( −r 2 r2 d ) exp ( −β (t−2tp tp )2 ) = IΩ(r,z)It(t) (2.7) where IΩ(r,z) = I0 1−Rf δ exp ( − r2 r2 d ) exp ( − z δ ) It(t)= √ β tp √ π exp ( −β (t−2tp tp )2 ) (2.8) 842 M. Ciesielski and I0 is laser intensity,Rf is reflectivity of the irradiated surface, δ is optical penetration depth, rd is characteristic radius of Gaussian laser beam, β=4ln2 and tp is characteristic time of the laser pulse. Here, it is assumed that the total time of the laser action beam on the surface is equal to 4tp. So, action of the laser beam is taken into account by introduction of the internal heat source Q(r,z,t). At the same time the dimensions Z and R are large enough that on the appropriate boundaries adiabatic conditions qb(r,z,t) = 0 can be assumed. In the case of the problem considered (see: Eq. (2.6)) one has (r,z)∈Γ : −λ [ n ·∇T(r,z,t)+ τT ∂[n ·∇T(r,z,t)] ∂t ] =0 (2.9) The initial conditions (the initial temperature of domain T0(r,z) and the initial heating rate v0(r,z) are also given t=0 : T(r,z,0)=T0(r,z) ∂T(r,z,t) ∂t ∣ ∣ ∣ ∣ ∣ t=0 = v0(r,z) (2.10) 3. Numerical solution using the Control Volume Method To solve the problem presented in the previous Section, the control volume method (CVM) is used. Thismethod constitutes a very effective tool for numerical modeling of heat transfer pro- cesses described by theFourier-type equations. In the case of numerical simulation ofmicroscale heat transfer and the models based on the DPL equation, this method has so far been applied only to the numerical solution using an ‘explicit’ scheme (Mochnacki and Ciesielski, 2015). The first stage of themethod application is the division of the domain considered into small cells (known as the control volumes CV). In this work, the shape of control volumes is regular one (it corresponds to the rings of a rectangular cross-section). Themore complex discretization using e.g. theVoronoi polygons can be also taken into account (Ciesielski andMochnacki, 2014). In Fig. 2, the domain discretization is presented, while in Fig. 3 the selected internal and boundary (top) control volumes are shown. Fig. 2. Discretization of the domain Application of the alternating direction implicit method... 843 Fig. 3. The internal and boundary control volumes On the basis of simple geometrical considerations, one can determine values of the successive volumes ∆Vi,j and surfaces (∆Ak)i,j limiting ∆Vi,j in each k-direction. Numerical modelling of transient problems requires introduction of the time grid, too: 0= t0 0 0 if j=0 (3.7) (q2)i,j =    (λ2)i,j [Ti+1,j −Ti,j ∆r + τT d dt (Ti+1,j −Ti,j ∆r )] if i 0 0 if i=0 (3.10) Application of the alternating direction implicit method... 845 and (λk)i,j are harmonic mean thermal conductivities between two central points of adjoining control volumes (λ1)i,j = 2λi,jλi,j−1 λi,j +λi,j−1 (λ2)i,j = 2λi,jλi+1,j λi,j +λi+1,j (λ3)i,j = 2λi,jλi,j+1 λi,j +λi,j+1 (λ4)i,j = 2λi,jλi−1,j λi,j +λi−1,j (3.11) and next, the thermal resistances are defined as follows (R1)i,j = ∆z (λ1)i,j (R2)i,j = ∆r (λ2)i,j (R3)i,j = ∆z (λ3)i,j (R4)i,j = ∆r (λ4)i,j (3.12) Then, Eq. (3.6) takes the form ∫ Ai,j [ n ·λ ( ∇T(r,z,t)+ τT ∂∇T(r,z,t) ∂t )] dA∼= 4 ∑ k=1 (θk)i,j (Rk)i,j (∆Ak)i,j (3.13) where (θ1)i,j = ( Ti,j−1−Ti,j + τT d(Ti,j−1−Ti,j) dt ) ∣ ∣ ∣ ∣ ∣ if j>0 (θ2)i,j = ( Ti+1,j −Ti,j + τT d(Ti+1,j −Ti,j) dt ) ∣ ∣ ∣ ∣ ∣ if i0 (3.14) while the notation expression|if condition introduced above, means expression ∣ ∣ ∣ if condition = { expression if condition = true 0 otherwise (3.15) After the introduction of all discrete terms into equation (3.1), one obtains ci,j (dTi,j dt + τq d2Ti,j dt2 ) ∆Vi,j = 4 ∑ k=1 (θk)i,j (Rk)i,j (∆Ak)i,j + ( Qi,j + τq dQi,j dt ) ∆Vi,j (3.16) or ci,j (dTi,j dt + τq d2Ti,j dt2 ) = 4 ∑ k=1 (θk)i,j (Rk)i,j (Φk)i,j +Qi,j + τq dQi,j dt (3.17) where (Φk)i,j =(∆Ak)i,j/∆Vi,j. 846 M. Ciesielski 3.2. Integration of the equation with respect to time The second stage of CVM is integration of equation (3.17) with respect to time. The same effect can be obtained introducing the approximation of time derivatives occurring in (3.17) by appropriate finite differences. The idea of the ADImethod is to split the time step∆t= tf+1− tf into two half-steps and apply two different finite difference schemes for each half time step. In the first half time step, a simple implicit scheme for directions (d1,d2) is used and simultaneously an explicit scheme for directions (d3,d4) is applied.Next, in the secondhalf time step, thedifference schemasarewritten by reversing the directions of the explicit and implicit schemes. The notation ‘(d1,d2)-(d3,d4)’, where the indexes di indicate directions of the neighbouring CV (see Fig. 1), is introduced. For passing: tf → tf+0.5 → tf+1, f =1, . . . ,F, and using the variant of ADI: (1,2)-(3,4), the following differential schemas are proposed c f i,j ( T f+0.5 i,j −T f i,j 0.5∆t + τq T f+0.5 i,j −2T f i,j +T f−0.5 i,j (0.5∆t)2 ) = ∑ k=1,2 (θk) f+0.5 i,j (Rk) f i,j (Φk)i,j + ∑ k=3,4 (θk) f i,j (Rk) f i,j (Φk)i,j +Q f+0.5 i,j +τq Q f+0.5 i,j −Q f i,j 0.5∆t (3.18) and c f+0.5 i,j ( T f+1 i,j −T f+0.5 i,j 0.5∆t + τq T f+1 i,j −2T f+0.5 i,j +T f i,j (0.5∆t)2 ) = ∑ k=1,2 (θk) f+0.5 i,j (Rk) f+0.5 i,j (Φk)i,j + ∑ k=3,4 (θk) f+1 i,j (Rk) f+0.5 i,j (Φk)i,j +Q f+1 i,j + τq Q f+1 i,j −Q f+0.5 i,j 0.5∆t (3.19) and (θk) s i,j for s∈{f,f+0.5,f +1} for this method are approximated in the following way (θ1) s i,j = [ Tsi,j−1−Tsi,j + τT ( Tsi,j−1−T s−0.5 i,j−1 0.5∆t − Tsi,j −T s−0.5 i,j 0.5∆t )]∣ ∣ ∣ ∣ ∣ if j>0 (θ2) s i,j = [ Tsi+1,j −Tsi,j + τT ( Tsi+1,j −T s−0.5 i+1,j 0.5∆t − Tsi,j −T s−0.5 i,j 0.5∆t )] ∣ ∣ ∣ ∣ ∣ if i0 (3.20) After transformations, the first system of equations (3.18) can be written in the final form as (A′0) f i,jT f+0.5 i,j +(A ′ 1) f i,jT f+0.5 i,j−1 ∣ ∣ ∣ if j>0 +(A′2) f i,jT f+0.5 i+1,j ∣ ∣ ∣ if i0 +(A′2) f i,j ∣ ∣ ∣ if i0 +2µT ( T f i,j −T f i+1,j )(Φ2)i,j (R2) f i,j ∣ ∣ ∣ ∣ ∣ if i0 +2c f i,j (1+4µq)T f i,j −2µqT f−0.5 i,j ∆t +(1+2µq)Q f+0.5 i,j −2µqQ f i,j (3.22) while the second system of equations (3.19) – in the following form (A′′0) f+0.5 i,j T f+1 i,j +(A ′′ 3) f+0.5 i,j T f+1 i,j+1 ∣ ∣ ∣ if j0 =(D′′) f+0.5 i,j (3.23) where (A′′k) f+0.5 i,j =−(1+2µT) (Φk)i,j (Rk) f+0.5 i,j k=3,4 (A′′0) f+0.5 i,j =2c f+0.5 i,j 1+2µq ∆t − ( (A′′3) f+0.5 i,j ∣ ∣ ∣ if j0 ) (D′′) f+0.5 i,j =2µT ( T f+0.5 i,j −T f+0.5 i,j+1 ) (Φ3)i,j (R3) f+0.5 i,j ∣ ∣ ∣ ∣ ∣ if j0 + [ −(1+2µT) ( T f+0.5 i,j −T f+0.5 i,j−1 ) +2µT ( T f i,j −T f i,j−1 ) ] (Φ1)i,j (R1) f+0.5 i,j ∣ ∣ ∣ ∣ ∣ if j>0 + [ −(1+2µT) ( T f+0.5 i,j −T f+0.5 i+1,j ) +2µT ( T f i,j −T f i+1,j ) ] (Φ2)i,j (R2) f+0.5 i,j ∣ ∣ ∣ ∣ ∣ if i