Transcript daycinmytalk - University of Kentucky

A Simple Numerical Approach For Solving A DualPhase-Lag Micro scale Heat Transport Equation
Illayathambi Kunadian
J. M. McDonough
Ravi Ranjan Kumar
Department of Mechanical Engineering
University of Kentucky, Lexington, KY 40506
Advanced Computational Fluid Dynamics
CF
1
D
Overview of this talk
I.
Introduction
II.
Brief review of origins of DPL equation
III. Discretization and analysis of unsplit 1-D DPL equation
I.
Stability analysis
IV.
Numerical scheme for solving 3-D DPL equation
V.
Computed results from selected problems
VI.
I.
1-D problem
II.
3-D problem
Summary and conclusions
Advanced Computational Fluid Dynamics
CF
2
D
Challenges in nanoscale heat transfer

From a microscopic point of view, ultrafast laser heating of metals is
composed of three processes:



deposition of radiation energy on electrons
transport of energy by electrons
heating of the material lattice through phonon-electron
interactions.

During a relatively slow heating process, the deposition of radiation energy
can be assumed to be instantaneous and can be modeled by Fourier
conduction; but applicability of this approach to very short-pulse laser
applications becomes questionable.

We must look for non-Fourier models because the laser pulse duration is
shorter than the thermalization time (time required for the phonons and
electrons to come into thermal equilibrium) and relaxation time of the
energy carriers.
Advanced Computational Fluid Dynamics
CF
3
D
Challenges in nanoscale heat transfer

An alternative is the hyperbolic heat conduction model; but this suffers
from violation of the second law of thermodynamics, and physically
unrealistic solutions are therefore unavoidable.

Successful attempts to model microscale heat transfer have been made by
Qui and Tien, but when investigating macroscopic effects a different
model is required.

Tzou proposed the dual phase lag model that reduces to parabolc,
hyperbolic, phonon-electron inteaction and pure phonon scattering models
under special values of relaxation times.
Advanced Computational Fluid Dynamics
CF
4
D
Origin of Dual phase lag model
 

q (r , t   q )  kT (r , t   T )
Tzou (1995)
T ~ delay behavior in establishing the temperature gradient
q ~ delay behavior in heat-flow departure
 

q (r , t )
[T (r , t )]
 


q (r , t )   q
 k T (r , t )   T

t
t


Energy equation
DPL
model

q
T
 S  C p
t
t
 q  2T 1 T
 ( 2T )
1
S 
2





T

S




T
q
 t 2  t
t
k
t 
Advanced Computational Fluid Dynamics
5
CF
D
Numerical Methods

Explicit Methods - 3D

Dai and Nassar developed implicit finite difference scheme


Split DPL equation into system of 2 equations and individual equations
solved using Crank-Nicolson scheme and solved sequentially

Discrete energy method to show unconditional stability of numerical
scheme
Zhang and Zhao employed iterative techniques like Gauss-Seidel, SOR, CG,
PCG to solve 3-D DPL equation

Used Dirichlet conditions, but applying Neumann boundary conditions
result in non-symmetric seven banded positive semi-definite matirces
not suitable for iterative methods like CG and PCG
Advanced Computational Fluid Dynamics
6
CF
D
Origin of Dual phase lag model

Present method






Formulation based on unsplit DPL equation
Stability shown using von Neumann stability analysis
Extend to 3D
Douglas –Gunn time splitting and delta-form Douglas Gunn time Splitting
Performance compared with numerical techniques available in literature
Results from specific problems
Advanced Computational Fluid Dynamics
7
CF
D
Laser heating source term
Gold Film
Laser
Tzou (1995)
Intensity of
L = 100 nm
laser absorption
1.992
1  R 
S 0  0.94 J 

t

 p 
I (t )  I 0 e

tp
(Intensity of laser )
S ( x, t )  S0e  I (t )
Qui and Tien (1992)
 t 
  2.77 
1  R  
tp 
 
S ( x, t )  0.94 J 
e
x
 t p 
x
t








J = 13.4 Jm2
R=0.93
t p=96fs
 = 15.3nm
 = 1.2104m2s1
q = 8.5ps
 = 90ps
T
k = 315Wm1K1
x 1.992 t  2t p
 
1  R  
tp
S ( x, t )  0.94 J 
e
 t p 
2
2
 t 
 2.77 
tp 
  (Intensity of laser
I (t )  I 0 e
2
 L x   L y
 x
  y
2  
2


1  R 
2ro2
S (r , t )  0.94 J 
e
)
2



  z  1.992 t  2t p

tp
 t p 
Advanced Computational Fluid Dynamics
3-D laser source
CF
8
D
Heat Conduction in a solid bar
Initial Conditions
 q  2T 1 T
 3T
 2T

T 2  2
2
 t
 t
x t x
T ( x,0)  T0
T
( x,0)  0
t
Boundary Conditions
T ( x, t )  T0
 ( x, t ) 
TW  T0

x
 q
Z 0
 2T T
 3T
 2T

Z 2  2
2
t
t
x t x
Z 1
 2T T    2T T 



0


2
2

t

t

t
x
 x


T (0, t )  TW
t
T
( x, t )  0
x
q
Z
t 0
T
q
x 
Initial Conditions
T ( x,0)  1 T ( x,0)  0
t
Boundary Conditions
T (0, t )  1
T
( x, t )  0
x
Advanced Computational Fluid Dynamics
t 0
x 
CF
9
D
Discretization and analysis of unsplit DPL equation
Trapezoidal integration
T  2T
 3T
 2T

Z

2
t t 2
x t 51×51×11
x2
n 1
n
 T 
 T 
Tmn 1  Tmn   
 
 t  m
 t  m
n 1
 T 
 
 t  m
Z
T
q
 2  n 1  2  n 
 2  n 1  2  n 
 T
 T   t   T 
 T 
 Z 



 x 2   2  x 2 
 x 2  
 x 2 







 



n
1
 T 
n 1
n 1  2T
1
1
n 1
n
n 1 

[
T

T
m
m1 ]

[3Tm  4Tm  Tm1 ] t 

[T
 2Tm  Tm 1 ]
  m 2t
2
2 m 1
2t
x
x

C4 Tmn11
 Tmn11

 C5Tmn 1

 C6 Tmn1
 Tmn1

 C7Tmn
1 n 1
 Tm
t
t  1 

C4   Z   2 
2  x 

t   1 

C6    Z    2 
2   x 

t   2  
1

C5   Z    2   1  
2   x   t 

t   2  
2

C 7    Z     2   1  
2   x   t 

Advanced Computational Fluid Dynamics
CF
10
D
Stability analysis




C4 Tmn11  Tmn11  C5Tmn 1  C6 Tmn1  Tmn1  C7Tmn 
1 n 1
Tm
t
Vmn 1  Tmn




C4 Tmn11  Tmn11  C5Tmn 1  C6 Tmn1  Tmn1  C7Tmn 
Tmn 1 
1 n
Vm
t
C6 2 cos h  C7 T n 
1
Vmn
m
C4 2 cos h  C5 
C4 2 cos h  C5 t
Tmn1  eihTmn
Tmn1  eihTmn
Vmn 1  Tmn
Tmn 1   C6 2 cos h   C7 

  C4 2 cos h   C5 


V n 1  
1
 m  

 

1
 T n 
m
C4 2 cos h   C5 t   
 
 n 
 Vm 
0
  

zmn 1  Cz mn
Advanced Computational Fluid Dynamics
CF
11
D
Stability analysis
Von Neumann
necessary condition for
Stability
 C6 2 cos h   C7 
 C 2 cos h   C 
5
 4
C

1


2 
 

1

C4 2 cos h   C5 t 


0


C 1
C6 2 cos h  C7   
1
0
C4 2 cos h  C5  C4 2 cos h  C5 t
C6 2 cos h   C7  
C4 2 cos h   C5 
 C6 2 cos h   C7  
1

  4
C4 2 cos h   C5 t
 C4 2 cos h   C5  
2
2
max   ,    1
Advanced Computational Fluid Dynamics
CF
12
D
Stability analysis
Distribution of

Advanced Computational Fluid Dynamics
CF
13
D
Finite difference scheme 3D
 q  2T 1 T
 ( 2T )
1
S 
2





T

S




T
q
 t 2  t
t 51×51×11 k 
t 
x 1.992 t  2t p
 
1  R  
tp
S ( x, t )  0.94 J 
e
t

 p 
2
 L x   L y
 x
  y
2  
2


1  R 
2ro2
S (r , t )  0.94 J 
e
2



  z  1.992 t  2t p

tp
 t p 
Advanced Computational Fluid Dynamics
CF
14
D
Finite difference scheme 3D
Trapezoidal integration
 q  2T 1 T
 ( 2T )
1
S 
2





T

S




T
q
 t 2  t
t 51×51×11 k 
t 
  2T  n 1   2T  n 
  2T  n 1   2T  n 




  2  
  2  
2 
 x 2 






 x 
 x 
 x  
 q n 1

 S in, j ,1k  S in, j , k

n




n

1
n
n

1
n
 Ti , j , k  Ti , j , k



2
  2T   t    2T 
  2T   t 
2

      T 

  2      2 
  2   
 1  T  n 1  T  n   T   y 2 
n

1
n
Si, j , k  Si, j , k 
 

 y   2   y 
 y   k 
 
  



q


n 1
n
n 1
n
   t  i , j , k  t  i , j , k  


2
2
2
2
2
 T  
 T  
  T 
  T 
  2  
  2  
  2 
  2 

z

z

z


 

 z  
 
 


n 1
 2T
1
 T 

[3Ti n, j,1k  4Ti n, j , k  Ti n, j,1k ]
 
 t  i , j , k 2t
n
1
 T 

[Ti n, j,1k  Ti n, j,1k ]
 
 t  i , j , k 2t
x 2
 2T
y 2
 2T
z 2



1
x 2
1
y 2
1
z 2
[Ti 1, j , k  2Ti , j , k  Ti 1, j , k ]
[Ti , j 1, k  2Ti , j , k  Ti , j 1, k ]
[Ti , j , k 1  2Ti , j , k  Ti , j , k 1 ]
Advanced Computational Fluid Dynamics
CF
15
D
Finite difference scheme 3D






C4Tin, j,1k  C5 Tin1,1j , k  Tin1,1j , k  C6 Tin, j11, k  Tin, j11, k  C7 Tin, j,1k 1  Tin, j,1k 1  F n






F n  C8Tin, j , k  C9 Tin1, j , k  Tin1, j , k  C10 Tin, j 1, k  Tin, j 1, k  C11 Tin, j , k 1  Tin, j , k 1 
 1 q 
  2C5  C6  C7 
C4   



t


t   1 

C5   T    2 
2   x 

t   1 

C6   T    2 
2   y 

t   1 

C7   T    2 
2   z 

 q n 1
Ti , j , k  tG *
t
t   1
1
1   1 2 q 


C8  2   T    2  2  2    
2   x
y
z    t 

t   1 

C9     T    2 
2   x 

t   1 

C10     T    2 
2   y 

t   1 

C11     T    2 
2   z 

 Sin, j ,1k  Sin, j , k
Sin, j ,1k  Sin, j , k 
G 
 q

2
t


*
Advanced Computational Fluid Dynamics
CF
16
D
Finite difference scheme 3D
 1 q 
 




t


1  2C
*
5







 2C6*  2C7* Tin, j,1k  C5* Tin1,1j , k  Tin1,1j , k  C6* Tin, j11, k  Tin, j11, k  C7* Tin, j,1k 1  Tin, j,1k 1  S n
S 
n
1  2C T
* n 1
5 i, j , k




 C5* Tin1,1j , k  Tin1,1j , k  S n
 2C6*Tin, j,1k  C6* Tin, j11, k  Tin, j11, k  0


 2C7*Tin, j,1k  C7* Tin, j,1k 1  Tin, j,1k 1  0
Advanced Computational Fluid Dynamics
Fn
 1 q 
 

  t 
C5
C5* 
 1 q 
 

  t 
C6
C6* 
 1 q 
 

  t 
C7
C7* 
 1 q 
 

  t 
CF
17
D
Finite difference scheme 3D
I  Ax T (1)  S n  AyT n  AzT n
I  Ax T (1)  S n  I  AT n
I  Ay T (2)  T (1)  AyT n
I  Ay T (2)  T (1)
I  Az T
I  Az T (3)  T (2)
T
Douglas-Gunn timesplitting
(3)
n 1
T
( 2)
T
 AzT
n
T n 1  T (3)  T n
(3)











I  Ax  1  2C5* Tin, j,1k  C5* Tin1,1j , k  Tin1,1j , k
I  Ay  1  2C6* Tin, j,1k  C6* Tin, j11, k  Tin, j11, k
I  Az  1  2C7* Tin, j,1k  C7* Tin, j,1k 1  Tin, j,1k 1
delta-form DouglasGunn time-splitting

Advanced Computational Fluid Dynamics
CF
18
D
Results
Z=0
Hyperbolic
Z = 0.01
Temp. Grad.
Precedence
Advanced Computational Fluid Dynamics
CF
19
D
Results
Z=1
Parabolic
Z = 100
Heat flux.
Precedence
Short pulse laser heating on thin metal film1-D
Advanced Computational Fluid Dynamics
CF
20
D
Results
Advanced Computational Fluid Dynamics
CF
21
D
3-D Schematic of femtosecond laser heating of gold film
200nm laser
beam
Work piece-Gold
250nm
500nm
100nm
250nm
500nm
500nm
500nm
500nm
3-D schematic of laser heating of gold film at different locations
Advanced Computational Fluid Dynamics
CF
22
D
Results
DPL
Parabolic
DPL DPL
Parabolic
Parabolic
Parabolic
Hyperbolic
Parabolic
DPL
DPL
Parabolic
Hyperbolic
Temperature distribution at top surface of gold film predicted by different models
Advanced Computational Fluid Dynamics
CF
23
D
Results
At t = 0.3 ps
DPL
Parabolic
DPL DPL
Parabolic
Parabolic
DPL
Parabolic
Hyperbolic
Hyperbolic
Hyperbolic
At t = 0.9 ps
DPL
Parabolic
Hyperbolic
Hyperbolic
Temperature distribution at top surface of gold film predicted by different models
Advanced Computational Fluid Dynamics
CF
24
D
Temperature distribution cont.
At t = 1.56 ps
DPL
Parabolic
DPL
Hyperbolic
Parabolic
At t = 2.23 ps
DPL
Parabolic
Hyperbolic
Parabolic
Temperature distribution at top surface of gold film predicted by different models
Advanced Computational Fluid Dynamics
CF
25
D
Performance comparison
Dai and Nassar
Total
N = 21
Explicit Method
4.88
Gauss-Siedel
13.46
Implicit Conjugate Gradient 12.12
Method D-G Time splitting
9.75
Delta D-G
9.15
Numerical Schemes
Present
Numerical Schemes
Explicit Method
Gauss-Siedel
Implicit Conjugate gradient
Method
D-G Time splitting
Delta D-G
CPU time in seconds
N = 41 N = 51 N = 101
147.62 450.26 7920.00
175.10 415.86 7800.00
110.50 233.96 2733.43
82.29 166.90 1792.36
75.22 153.30 1637.96
Total CPU time in seconds
N = 21 N = 41 N = 51 N = 101
4.88 147.62
450.26 7920.00
14.14 253.42
627.03 11343.06
12.33 124.83
270.30 3614.69
9.24 82.44
165.76 1506.38
8.54 70.50
140.92 1344.40
Advanced Computational Fluid Dynamics
CF
26
D
Summary and conclusions
New numerical technique to solve DPL implicitly

Unconditionally stable numerical scheme for solving 1-D DPL equation


Solves one equation instead of splitting DPL equation into 2
equations and apply discretization

Reduces number of arithmetic operations involved

Reduces computational time
New formulation satisfies von Neumann necessary condition for
stability

Heat conduction in a solid bar

Semi-infinite slab – temperature raised at one end

q
is responsible for presence of sharp wave front in heat propagation in
CHE conduction

T
diminishes the sharp wave front and extends heat affected zone
deeper into the medium
Advanced Computational Fluid Dynamics
CF
27
D
Summary and conclusions
Numerical scheme for solving 3-D DPL equation


The new numerical formulation of discretizing DPL directly
outperforms Dai’s method of splitting DPL into two equations and
then apply discretization

Delta-form Douglas-Gunn time-splitting method outperforms all
other numerical techniques – CPU time taken for entire simulation

Explicit method good for small N (N=21). N > 21 all implicit
methods except Gauss-Seidel method perform better than
explicit method.

CV wave and diffusion models predict higher temperature level in heat
affected zone than the DPL model, but penetration depth is much
shorter - formation of thermally undisturbed zone.

DPL model - Heat affected zone is significantly larger than other models

Also, DPL results in 3D exhibit similar behavior as the one-dimensional
results
Advanced Computational Fluid Dynamics
CF
28
D