Modeling of laser heating GaAs considering the effects of atmospheric thermal blooming with crosswind

Modeling of laser heating GaAs considering the effects of atmospheric thermal blooming with crosswind

Optik 131 (2017) 11–20 Contents lists available at ScienceDirect Optik journal homepage: www.elsevier.de/ijleo Original research article Modeling ...

828KB Sizes 0 Downloads 21 Views

Optik 131 (2017) 11–20

Contents lists available at ScienceDirect

Optik journal homepage: www.elsevier.de/ijleo

Original research article

Modeling of laser heating GaAs considering the effects of atmospheric thermal blooming with crosswind Guibo Chen ∗ , Chaoqun Wang School of Science, Changchun University of Science and Technology, No.7089 of WeiXing Road, ChaoYang District, ChangChun City, JiLin Province, 130022, PR China

a r t i c l e

i n f o

Article history: Received 21 September 2016 Accepted 12 November 2016 Keywords: Laser heating GaAs Thermal blooming 2-D modeling

a b s t r a c t In this paper, a two-dimensional (2-D) modeling for temperature of GaAs irradiated by laser considering the effects of atmospheric thermal blooming with crosswind is presented. By using analytical expression of distorted laser intensity for the atmosphere thermal blooming with crosswind, the theoretical 2-D plane heat transfer model of laser irradiating GaAs after the action of atmosphere thermal blooming is established. Integral transform method is used to solve the heat transfer equation and its analytical solutions are obtained. The effects of thermal blooming on temperature distributions of GaAs are studied. Results show that the crosswind and the laser propagation distance have obvious influences on the thermal blooming, thereby affecting the temperature distributions of GaAs induced by distorted laser. The smaller the crosswind velocity is, the stronger the thermal blooming is, which leads to a higher distortion of temperature distributions. On the other hand, the farther the laser propagation distance in atmosphere is, the stronger the thermal blooming is, which also leads to the higher distortion of temperature distributions. © 2016 Elsevier GmbH. All rights reserved.

1. Introduction Modeling of laser heating is very important for laser applications such as laser processing and laser damage, because it can not only reduce the experimental cost and time, but also can provide results under many conditions including the extreme conditions that the experimentation is difficult to achieve. Moreover, analytical modeling has always been very attractive to researchers because it can establish an intuitive relationship between parameters and results, which has important significance to reveal the physical mechanism of laser-matter interaction and optimize the laser or workpiece parameters [1–8]. In the current studies of laser heating, most of laser beams are assumed to irradiate onto material surface directly without any medium. But in the potential applications of laser in the fields such as industry or military, it often encounters that the laser propagates in the atmosphere, and then irradiates onto the material surface. Laser propagation in the atmosphere is a very complicated physical process, especially the thermal blooming induced by laser heating air. When the laser propagates in the atmosphere, the nonlinear effects will cause the distortion of the laser beam, which will lead to the thermal blooming [9–15]. Thermal blooming can affect the quality of the laser beam, which leads to the instability of the laser propagation,

∗ Corresponding author. E-mail address: [email protected] (G. Chen). http://dx.doi.org/10.1016/j.ijleo.2016.11.078 0030-4026/© 2016 Elsevier GmbH. All rights reserved.

12

G. Chen, C. Wang / Optik 131 (2017) 11–20

and affects the interaction between the laser and material. Thermal blooming is very important for the study of intense laser propagation in the atmosphere, which dominates laser distortion due to atmosphere. After the effects of the atmospheric thermal blooming, the laser intensity distribution on the material surface is no longer regular. At this time, it needs to solve a 2-D or 3-D heat transfer problem when the irregular laser heating material is modeled. In this paper, we present 2-D plane modeling for distorted laser heating material considering the atmospheric thermal blooming with crosswind and investigate the effects of thermal blooming on the heat transfer in the material. GaAs is widely used in the optoelectronic devices because of its wide band gap and direct band gap. When the optoelectronic devices work under the intense laser environment, they often face performance degradation and failure even induced damage due to laser heating. It is important to study the interaction of laser with GaAs and its heating mechanisms for destruction and protection of the optoelectronic devices. So we present a 2-D modeling of laser heating GaAs considering the effects of atmospheric thermal blooming in this paper. It is organized as follows. In Section 2 we introduce the theoretical models of thermal blooming and distorted laser heating GaAs, and analytical solutions of the temperature distributions induced by a thermal blooming laser heating are obtained. Results of thermal blooming and temperature for different crosswind and propagation distance are presented in Section 3. Our main conclusions are summarized in Section 4. It is pointed out that, in this paper, a collimated continuous wave (CW) laser is considered, and the thermal blooming is assumed to be steady state. 2. Mathematical modeling 2.1. Thermal blooming modeling We assume that d is the laser propagation distance in the atmosphere and x,y are the variables in transverse plane. From reference [9,10], the scalar wave equation in parabolic approximation is,



∇ 2T

+ 2ik0



∂ n2 (x, y, d, t) −1 + k02 ∂z n20

=0

(1)

where, n is refractive index of atmosphere, n0 is static refractive index, k0 is wave number,

∇ 2T =

2

is the field amplitude and

2

∂ + ∂ . ∂x2 ∂ y2

The structural equation can be written as: n2 n20



 1

− 1 ≈ 2 n0 − 1

(2)

0

where, 0 is the density of the atmosphere with no distorting, and 1 is the change of the density of atmosphere and yields,

∂1 1− ˛Ip = ∂t cs2 where,  =

cp cv ,

(3)

cp and cv is specific heat capacity at constant pressure and constant volume respectively, cs2 =



 ∂p is the ∂

sound velocity, p is the atmospheric pressure, and p is the density of atmosphere. When the laser is a collimated beam, the diffraction effect can be neglected due to the large Rayleigh distance. Form the reference [16], in the steady state, the expression of the thermal blooming laser intensity can be written as:

 





Ip (x, y, d) = I0 e−d exp − x2 + y2 /r02 exp [−NC 0 (x, y) g 1 (d)]

(4)

where,  is the laser absorption in atmosphere, r0 is the waist radius of Gaussian beam, and Nc , 0 (x, y) , g 1 (d) yield: NC = −

dn d2 I0 dT cp vx r0

(5)

    2x exp − x2 + y2 /r02 + 0 (x, y) = r0 g1 (d) =

2 d

1−

1 − e−d d



1−

4y2 r02



y2

 − r2 e 0 2



1 + erf

x  r0

(6)

(7)

dn is the temperature change where, vx is the crosswind velocity along the x direction, NC is the distortion parameter, and dT ratio of the refractive index of the atmosphere. The distortion parameter NC is the key parameter to describe the intensity of thermal blooming, we will discuss the effect of distortion parameter on thermal blooming in the third part of this paper.

G. Chen, C. Wang / Optik 131 (2017) 11–20

13

Fig. 1. Schematic diagram of laser irradiating on GaAs after propagation through atmosphere.

2.2. Laser heating modeling If the thickness of the GaAs is very thin, temperature distributions in the thickness direction tend to be uniform quickly when the laser irradiating. At this time, the two-dimensional (2-D) plane heat transfer model can be established. As shown in Fig. 1, we assume that the size of GaAs bulk in two directions x, y in a 2-D Cartesian coordinate system are Lx , Ly respectively. Considering the body absorption mechanism of GaAs to the laser energy, 2-D classical Fourier heat conduction equation can be written as: 2

2

∂ T (x, y; t) ∂ T (x, y; t) Q (x, y; t) 1 ∂T (x, y; t) = + + ˛ k ∂x2 ∂ y2 ∂t

(8)

where, ˛ = k/c is the thermal diffusivity, k is the thermal conductivity of the GaAs,  is the mass density, c is the heat capacity, and Q (x, y; t) is heat source for GaAs absorption of laser energy,





Q (x, y; t) = 1 − rf ıIp (x, y)  (t)

(9)

where, rf is the reflection coefficient, ı is absorption coefficient of GaAs,  (t) is time distribution function of laser intensity and it yields  (t) = 1, Ip (x, y) is the laser intensity distribution function after thermal blooming as Eq. (4). It is assumed that the boundary conditions of Eq. (8) are adiabatic as:

∂T (x, y; t) ∂T (x, y; t) = =0 | | ∂x ∂y x=0,Lx x=0,Ly

(10)

and the initial condition yields: T (x, y; t) |t=0 = T0

(11)

In order to solve the Eqs. (8)–(11), forward and inverse transforms for T (x, y; t) about x, y are introduced based on integral transformation method [17] respectively and written as:

  ⎧ ∞    X ˇm , x ⎪ ⎪   · T¯ ˇm , y; t ⎪ ⎨ T (x, y; t) = m=1

 ⎪   ⎪ ⎪ ¯ ⎩ T ˇm , y; t =

N ˇm



Lx

X ˇm , x

x =0



 



(12)



T x , y; t dx



⎧ ∞     Y (n , y) ¯  ⎪ ⎪ ¯ ˇm , y; t = T · T¯ ˇm , n ; t ⎪ ⎨ n=1

⎪   ⎪ ⎪ ⎩ T¯¯ ˇm , n ; t = 





N (n )

Ly

y =0



Y n , y



 

(13)



T¯ ˇm , y ; t dy

where, X ˇm , x , Y (n , y) and ˇm ,  n are eigenfunctions and eigenvalues for the eigen problem obtained by separating





variables of the homogeneous part of the heat conduction equation (8) respectively, N ˇm , N (n ) are the normaliza-





tion integrals of ˇm ,  n respectively. We note that the eigenfunctions X ˇm , x , Y (n , y) of this eigen problem satisfy the orthogonal condition. T¯ and T¯¯ denote the transform with respect to the x and (x,y) respectively.

14

G. Chen, C. Wang / Optik 131 (2017) 11–20

Above problem can be solved by the successive application of the 1-D integral transform to the x and y variables, solving the resulting ordinary differential equation, and then inverting the transform of temperature successively. Then the solution of Eq. (8) can be written as:

  ∞ ∞   X ˇm , x Y (n , y)   2     T (x, y; t) = + n2 t × · exp −˛ ˇm 



N ˇm N (n )

m=1 n=1

 

t

T¯¯ 0 +

t  =0

  ˛   2 exp ˛ ˇm + n2 t  · Q¯¯ ˇm , n ; t  dt 



(14)

k

where, the double transforms T¯¯ 0 and Q¯¯ are defined as:



Lx

T¯¯ 0 =



0

Ly



0



 



X ˇm , x Y n , y T0 dx dy



Q¯¯ ˇm , n ; t  =



Lx



0

Ly



 

(15)

 



X ˇm , x Y n , y Q x , y ; t dx dy

(16)

0

According to the boundary condition, the eigenfunctions are:







X ˇm , x = cos ˇm x



(17)

Y ( n , y) = cos ( n y) , and the normalization integrals:







N ˇm =

 N (n ) =

Lx forˇm = 0

(18)

Lx /2forˇm = / 0 Ly forn = 0

(19)

Ly /2forn = / 0









The eigenvalues ˇm ,n ,are the non-negative roots of equations sin ˇm Lx = 0 and sin n Ly = 0. Substituting the eigenfunctions (17)–(19) into Eq. (15), analytical expression of T¯¯ can be obtained:



T¯¯ 0 = T0 ·

sin ˇm Lx





·

ˇm

sin n Ly

0



(20)

n





Substituting the eigenfunctions (17)–(19) and expression of Q (x, y; t) as (9) into Eq. (16), the integral Q¯¯ ˇm , n ; t  can be written as:









Q¯¯ ˇm , n ; t  = 1 − rf ı˝1 where,



Lx



˝1 = 0

Ly





(21)



 



cos ˇm x cos n y I p x , y dx dy

(22)

0

We note that the double integral ˝1 needs to be calculated by numerical double integration, the Gauss–Kronrod quadrature method is used to handle it [18]. Substituting Eqs. (15)–(22) into Eq. (14), The temperature T (x, y; t) of any location (x, y) in the material can be calculated at any time t. 3. Results and discussions Temperature of GaAs irradiated by CW laser considering the effects of atmospheric thermal blooming is calculated using above 2-D analytical solutions. The thermo-physical and geometric parameters of GaAs are given in Table 1, and the atmospheric parameters are shown in Table 2. We assume that the laser wavelength is 1.06mm, the laser peak power intensity is 106 W/m2 , the laser waist radius is 2 cm, and the initial temperature of GaAs is 300 K. All of the following calculated examples are using the above computational parameters. 3.1. Method validation It is assumed that the collimation laser is irradiated to the center of the GaAs surface (as shown in Fig. 1), the propagation distance of laser in atmosphere is 1.2 km, the crosswind velocity along x direction is 2m/s, and the laser irradiating time is 2 ms. The temperature distributions for fixed y = Ly /2 and different x were calculated by using the analytical solutions (14) of

G. Chen, C. Wang / Optik 131 (2017) 11–20

15

Table 1 Thermo-physical and geometric parameters of GaAs. Parameters

Values

density/(kg/m3 ) specific heat/(J/kg K) coefficient of heat conductivity/(W/m K) reflection coefficient absorption coefficient/m length of material/cm width of material/cm

5320 318 42.5 0.5 106 10 10

Table 2 Atmospheric parameters. Parameters

Values

absorption coefficient/m temperature change ratio of the refractive index density/(kg/m3 ) specific heat/(J/kg K)

10−5 −10−7 1.34 1000

800

T/K

600

M=N=5 M=N=10 M=N=20 FEM

400

200 2

3

4

5

6

7

x/cm Fig. 2. Temperature distributions calculated by analytical solutions (different M,N) and FEM.

this paper and the finite element method (FEM) of the published literature [2] as shown in Fig. 2. We note that in expression (14), only the finite M,N terms in the series of m,n need to be evaluated, while the rest of the terms can be omitted. It can be seen that the analytical solutions of this paper converges quickly with the gradual increase of M,N, and the calculated results of the analytical solutions and the FEM are very close when M,N are equal to 20. Therefore, we can conclude that the presented method is correct and effective. Taking into account the calculation accuracy and computational efficiency, the following examples in this paper are calculated based on M = N = 30. In addition, we have also done a lot of validation samples, and results show that the analytical method is correct. Due to limited space, we cannot present all the results. 3.2. Results and discussions We assume that the propagation distance of the laser in the atmosphere is 1.5 km, and the laser irradiating time is 2 ms. Figs. 3, 5, 7 and 9 give the laser intensity distributions when the crosswind velocity is 0.75, 1.5, 3 and 5 m/s respectively. It can be seen that the laser beam will be distorted and form a “crescent-shaped” distribution bending to headwind direction when there is a certain crosswind. This is because that air density decreases in the optical axis along the direction of the wind, so the refractive index decreases in the direction of the headwind; it causes the distortion and bending of the beam, and affects the quality of the beam and the laser power density at the target surface. In addition, it can be seen that crosswind velocity has a significant effect on the thermal blooming of laser propagation in the atmosphere. The larger the crosswind velocity is, the smaller the thermal blooming effect is; and on the contrary, the smaller the crosswind velocity is, the stronger the thermal blooming effect is. This is because the large crosswind velocity can quickly take the heat away from the optical path, so that the distortion of refractive index of the optical path becomes smaller. We can also explain this by distortion parameter NC as (5). Distortion parameters NC for vx = 0.75, 1.5, 3, 5 m/s are 1.12, 0.56, 0.28, 0.16 respectively, and the larger the distortion parameter, the stronger the thermal blooming effect is. Therefore, the crosswind velocity is helpful to reduce the thermal blooming effect of laser propagation in the atmosphere.

16

G. Chen, C. Wang / Optik 131 (2017) 11–20 10 1.2E6

8

1.05E6 9E5 7.5E5

y/cm

6

6E5 4.5E5

4

3E5 1.5E5

2

0

0

0

2

4

6

8

10

x/cm Fig. 3. Laser intensity distributions on the material surface when crosswind velocity is 0.75m/s. 10 900.0

8

825.0 750.0 675.0

y/cm

6

600.0 525.0

4

450.0 375.0

2

0

300.0

0

2

4

6

8

10

x/cm Fig. 4. Temperature distributions on the material surface when crosswind velocity is 0.75m/s. 10 1E6

8

8.75E5 7.5E5 6.25E5

y/cm

6

5E5 3.75E5

4

2.5E5 1.25E5

2

0

0

0

2

4

6

8

10

x/cm Fig. 5. Laser intensity distributions on the material surface when crosswind velocity is 1.5m/s.

Next, we calculate the temperature rise of GaAs induced by distorted laser radiation after propagation through the atmosphere, as shown in Figs. 4, 6, 8 and 10 due to the crosswind velocity of 0.75, 1.5, 3 and 5 m/s respectively. It can be seen that the shape of the temperature distributions are similar to that of the distorted laser by contrast profiles of laser intensity and temperature. This is due to its energy gain mechanism and the temperature rise is mainly dependent on the direct absorption of laser energy on the material surface. In addition, it finds that when the crosswind velocity is small, the beam distortion is larger, which leads to the larger temperature distortion. Next, we assume that the crosswind velocity is 1m/s, and the laser irradiating time is 2 ms. The laser intensity and the GaAs temperature distributions are calculated when the laser propagation distances are 1, 1.4 and 1.8 km respectively, as shown in Figs. 11–16 . It can be seen that the laser propagation distance has a significant effect on the thermal blooming. The larger the laser propagation distance is, the stronger the thermal blooming effect is. This is due to the larger laser propagation distance, the greater the negative lens effect of the atmosphere, and it leads to the bigger beam distortion. At the same time, the larger laser propagation distance, the more the thermal distortion accumulation of the beam, this also leads to the bigger beam distortion. We can also explain this by distortion parameter NC as (5). Distortion parameters NC for d = 1, 1.4, 1.8 km are 0.37, 0.73, 1.21 respectively. And the larger the distortion parameter, the stronger the thermal blooming effect is.

G. Chen, C. Wang / Optik 131 (2017) 11–20

17

10 800.0

8

737.5 675.0 612.5

y/cm

6

550.0 487.5

4

425.0 362.5

2

0

300.0

0

2

4

6

8

10

x/cm Fig. 6. Temperature distributions on the material surface when crosswind velocity is 1.5m/s. 10 9E5

8

7.875E5 6.75E5 5.625E5

y/cm

6

4.5E5 3.375E5

4

2.25E5 1.125E5

2

0

0

0

2

4

6

8

10

x/cm Fig. 7. Laser intensity distributions on the material surface when crosswind velocity is 3m/s. 10 750.0

8

693.8 637.5 581.3

y/cm

6

525.0 468.8

4

412.5 356.3

2

0

300.0

0

2

4

6

8

10

x/cm Fig. 8. Temperature distributions on the material surface when crosswind velocity is 3m/s. 10 9E5

8

7.875E5 6.75E5 5.625E5

y/cm

6

4.5E5 3.375E5

4

2.25E5 1.125E5

2

0

0

0

2

4

6

8

10

x/cm Fig. 9. Laser intensity distributions on the material surface when crosswind velocity is 5m/s.

18

G. Chen, C. Wang / Optik 131 (2017) 11–20 10 800.0

8

737.5 675.0 612.5

y/cm

6

550.0 487.5

4

425.0 362.5

2

0

300.0

0

2

4

6

8

10

x/cm Fig. 10. Temperature distributions on the material surface when crosswind velocity is 5m/s. 10 1E6

8

8.75E5 7.5E5 6.25E5

y/cm

6

5E5 3.75E5

4

2.5E5 1.25E5

2

0

0

0

2

4

6

8

10

x/cm Fig. 11. Laser intensity distributions on the material surface when propagation distance is 1 km. 10 800.0

8

737.5 675.0 612.5

y/cm

6

550.0 487.5

4

425.0 362.5

2

0

300.0

0

2

4

6

8

10

x/cm Fig. 12. Temperature distributions on the material surface when propagation distance is 1 km. 10 1E6

8

8.75E5 7.5E5 6.25E5

y/cm

6

5E5 3.75E5

4

2.5E5 1.25E5

2

0

0

0

2

4

6

8

10

x/cm Fig. 13. Laser intensity distributions on the material surface when propagation distance is 1.4 km.

G. Chen, C. Wang / Optik 131 (2017) 11–20

19

10 800.0

8

737.5 675.0 612.5

y/cm

6

550.0 487.5

4

425.0 362.5

2

0

300.0

0

2

4

6

8

10

x/cm Fig. 14. Temperature distributions on the material surface when propagation distance is 1.4 km.

10 1.4E6

8

1.225E6 1.05E6 8.75E5

y/cm

6

7E5 5.25E5

4

3.5E5 1.75E5

2

0

0

0

2

4

6

8

10

x/cm Fig. 15. Laser intensity distributions on the material surface when propagation distance is 1.8 km.

10 1000

8

912.5 825.0 737.5

y/cm

6

650.0 562.5

4

475.0 387.5

2

0

300.0

0

2

4

6

8

10

x/cm Fig. 16. Temperature distributions on the material surface when propagation distance is 1.8 km.

In addition, it can be seen that the temperature of GaAs also has some distortion, and the shape of temperature distributions is similar to that of the laser intensity distributions. The greater the laser propagation distance in the atmosphere, the stronger the thermal blooming effect, which leads to the greater temperature distortion of the GaAs.

4. Conclusions (1) Analytical solutions of 2-D modeling for temperature of GaAs irradiated by distorted laser considering the effects of atmospheric thermal blooming with crosswind are obtained by using 2-D integral transform method. It can be used to solve the 2-D plane heat transfer problem induced by the distortion source. (2) The crosswind has an obvious influence on the thermal blooming, thereby affecting the temperature distribution of GaAs. The smaller the crosswind velocity is, the stronger the thermal blooming is, which leads to a higher distortion of temperature of GaAs.

20

G. Chen, C. Wang / Optik 131 (2017) 11–20

(3) The laser propagation distance has an obvious influence on the thermal blooming, thereby affecting the temperature distribution of GaAs. The larger the laser propagation distance is, the stronger the thermal blooming is, which leads to a higher distortion of temperature of GaAs. Acknowledgement This work was supported by the Changchun University of Science and Technology Young Scientists Project “Numerical modeling and inverse problems of heat transfer in laser induced materials” under Grant No. XQNJJ-2014-03. References [1] B.S. Yilbas, Laser Heating Applications: Analytical Modeling, Elsevier, Amsterdam, 2012. [2] J. Bi, X.H. Zhang, X.W. Ni, G.Y. Jin, C.L. Li, Numerical simulation of thermal damage effects on Gallium Arsenide (GaAs) induced by a 0.53 ␮m wavelength long pulsed laser, Lasers Eng. 22 (2011) 37–46. [3] J. Bi, G.B. Chen, Characteristic analysis of 532 nm millisecond pulse laser-induced surface damage in the shape of thermal decomposition to GaAs by means of a Semi-analytical method, Lasers Eng. 32 (2015) 99–117. [4] Y. Qin, Y.B. Chen, X.W. Ni, Z.H. Shen, J. Bi, X.H. Zhang, Axisymmetric numerical simulation of plastic damage in aluminum alloy induced by long pulsed laser, Opt. Laser Eng. 48 (2010) 361–367. [5] R. Gospavic, M. Sreckovic, V. Popov, Modelling of laser-material interaction using semi-analytical approach, Math. Comput. Simulat. 65 (2004) 211–219. [6] Z.W. Li, H.C. Zhang, Z.H. Shen, X.W. Ni, Time-resolved temperature measurement and numerical simulation of millisecond laser irradiated silicon, J. App. Phys. 114 (2013) 033104. [7] M. Kalyon, B.S. Yilbas, Repetitive laser pulse heating analysis: pulse parameter variation effects on closed form solution, Appl. Surf. Sci. 252 (2006) 2242–2250. [8] B.S. Yilbas, A closed form solution for temperature rise inside solid substrate due to time exponentially varying pulse, Int. J. Heat Mass Trans. 45 (2002) 1993–2000. [9] J.A. Fleck, J.R. Morris, M.D. Feit, Time-dependent propagation of high-energy laser beams through the atmosphere: II, Appl. Phys. A 14 (1978) 99–115. [10] J.A. Fleck, J.R. Morris, M.D. Feit, Time-dependent propagation of high energy laser beams through the atmosphere, Appl. Phys. 10 (1976) 129–160. [11] G. Frederick, Twenty-five years of thermal blooming:an overview, SPIE 1221 (1990) 2–25. [12] Julians Nichols, Thermal blooming induced low order aberrations, Appl. Opt. 23 (1984) 2342–2347. [13] S. Enguehard, Review of the physics of small scale thermal blooming in up link propagation, SPIE 1991 (1415) 128–137. [14] V.F. Boris, P.L. Vladimir, Estimation of turbulent and thermal blooming degradation and required characterization of adaptive system, SPIE 3706 (1999) 361–367. [15] Yinbo Huang, Yingjian Wang, The scaling laws of laser beam spreading induced by turbulent and thermal blooming, SPIE 5639 (2004) 65–69. [16] P. Hillion, S. Quinnez, Thermal blooming calculations with analytical diffraction approximated expressions, J. Math. Phys. 22 (1981) 897–907. [17] M.N. Özisik, Heat Conduction, John Wiley and Sons, New York, 1993. [18] D.P. Laurie, Calculation of Gauss-Kronrod quadrature rules, Math. Comput. 66 (1997) 1133–1145.