Analytical model to investigate distributions of the thermal stresses in the pad and disk for different temporal profiles of friction power

Introduction

Phenomenon which determines friction, wear and reliability of frictional elements of the machines working in the braking mode occurs due to high temperature and also significant thermal stresses. Thermal stresses are caused by simultaneous effects of the temperature gradients and their maximum values^1,2 which depend on factors such as a braking mode, microgeometry, and macrogeometry of the contact surfaces and physicochemical properties of the materials. ³

In practice, there are known cases, where only the gradients of temperature cause so high level of stresses in the subsurface layers of a pad and a disk, that the friction process takes place in the plasticity contact conditions, even if the volumetric temperature of the friction elements and the contact pressure is so low, that in normal conditions, it could result in the elastic deformations only.^4,5 High values of the thermal stresses are conducive to the cracking initiation that gradually leads to disruption of the friction node work.^6–8 Thus, the estimation of the stresses distribution in the frictional couple elements of a braking system is important. The knowledge about the state of stresses allows to point the manner of stress value reduction and thus decrease the abrasive wear of friction surface and increase the expected exploitation time of the friction node.

Based on analysis of the results presented in the literature,^9–11 the topic of mathematical modeling and experimental investigations of the distributions of temperature and thermal stresses in the friction elements of a braking system is an intensively developing scientific problem. Theoretical works related to this problem can be divided into two main groups depending on applied methods of solution to the thermal problems of friction—the boundary-value problems of heat conduction taking into account the heat generation due to friction on the contact surface of the friction couple elements and the corresponding boundary problems of thermoelasticity. The analytical and analytical–numerical methods belong to the first group, and numerical methods belong to the second one.^12,13

Issues related to determination of the quasi-static thermal stresses in the friction elements of a braking system, based on exact solutions to the one-dimensional thermal problems of friction, were considered in the literature.^14–18 Calculation schemes devised in those works concern modeling of the frictional heating of the pad–disk system by means of appropriate solutions for two replacement systems: (1) two-element (strip–semi-space) and (2) three-element (consisting of a homogeneous half-space and a protective strip deposited on a semi-infinite homogeneous foundation). In both schemes, the case of a single braking process with a constant retardation was considered when the specific friction power linearly decreases with time from maximum value at the initial moment to zero at the moment of standstill. As experimental researches show, the temporal profile of specific friction power during braking can significantly differ from the linear one. Functions that describe the evolution of the specific friction power during braking were classified in the monograph ¹⁹ containing 10 different functions. Taking account of the profiles, at first, the influence of the specific friction power on the temperature during separate heating of one element of the friction couple was investigated. ²⁰ Next, the thermal problems of friction were considered. Temperature fields, obtained from the solution of the heat conduction problem for disk and pad (two semi-spaces), taking into account 3 of the 10 above-mentioned temporal profiles of specific friction power, were analyzed in Topczewska, ²¹ and the corresponding distributions of the thermal stresses were studied in Topczewska. ²² Exact equations to carry out the analysis of the temperature distributions in the pad and disk for all 10 time profiles of specific friction power were designated in Yevtushenko et al. ²³ Based on the obtained results, corresponding formulas to analyze the quasi-static thermal stresses were established in this article.

Temperature

Based on the tests on the thermal stability of the friction couples with simultaneous registration of time-varying braking torque, the following temporal profiles of the specific friction power were determined ²⁰

q i ( t ) = q 0 q i * ( t ) , q 0 = w 0 / t s , 0 ≤ t ≤ t s , i = 1 , 10 ¯

(1)

where t is the time, t s is the braking time, w 0 is the total friction work of the system at the stop time moment t = t s , which is equal for all 10 time profiles (1)

w 0 = ∫ 0 t s q i ( t ) dt , i = 1 , 10 ¯

(2)

and forms of dimensionless functions q i * ( t ) with the maximum values q i , max * , and the time of its achieving t i , max were presented in Table 1, where t i , max * = t i , max / t s and t * = t / t s .

OPEN IN VIEWER

Table 1.

Temporal profiles of specific friction power q i * ( t ) , i = 1 , 10 ¯ and their maximum values q i , max * achieved at the moments of time t i , max * .

i	q i * ( t )	q i , max *	t i , max *
1	2 ( 1 − t * )	2	0
2	2 t *	2	1
3	3 2 1 − t *	1.5	0
4	3 2 t *	1.5	1
5	3 ( 1 − t * ) 2	3	0
6	3 t * 2	3	1
7	6 t * ( 1 − t * )	1.5	0.5
8	6 5 ( 1 − t * ) ( 1 + 2 t * )	1.35	0.25
9	18 5 t * ( 1 − 2 3 t * )	1.35	0.75
10	6 t * ( 1 − t * )	1.5	0.25

To describe the frictional heating process in the pad–disk tribosystem during single braking, we choose the one-dimensional thermal problem of friction by taking into account the heat generation due to friction process on the contact surface of two different semi-spaces z ≥ 0 and z ≤ 0 (Figure 1). We assume that the initial temperature of both elements is the same and equal to the ambient temperature T a , and thermal contact on the sliding surface z = 0 is perfect

T 1 , i ( 0 , t ) = T 2 , i ( 0 , t ) , 0 < t ≤ t s

(3)

K 2 ∂ T 2 , i ( z , t ) ∂ z | z = 0 − K 1 ∂ T 1 , i ( z , t ) ∂ z | z = 0 = q i ( t ) , 0 < t ≤ t s

(4)

where the functions of specific friction power q i ( t ) , i = 1 , 10 ¯ have forms (1) and (2); K l are the coefficients of thermal conductivity of the pad ( l = 1 ) and disk ( l = 2 ) materials; and T l , i ( z , t ) is the transient temperature fields which disappear when the distance from the contact surface increases.

OPEN IN VIEWER

Figure 1.

Scheme of the problem.

Solutions to such formulated problem have the form^21,23

T l , i ( z , t ) = T a + T 0 T l , i * ( ζ , τ ) , 0 ≤ t ≤ t s ( 0 ≤ τ ≤ τ s ) , l = 1 , 2 , i = 1 , 10 ¯

(5)

where

T l , 1 * ( ζ , τ ) = 4 3 γ τ τ * 〈 { 3 / τ * − 2 [ 1 + Z l 2 ( ζ , τ ) ] } ierfc Z l ( ζ , τ ) + Z l ( ζ , τ ) erfc Z l ( ζ , τ ) 〉

(6)

T l , 2 * ( ζ , τ ) = 4 3 γ τ τ * { 2 [ 1 + Z l 2 ( ζ , τ ) ] ierfc Z l ( ζ , τ ) − Z l ( ζ , τ ) erfc Z l ( ζ , τ ) }

(7)

T l , 3 * ( ζ , τ ) = γ 〈 3 τ ierfc Z l ( ζ , τ ) − τ τ * { [ 1 + Z l 2 ( ζ , τ ) ] ierfc Z l ( ζ , τ ) − 0 . 5 Z l ( ζ , τ ) erfc Z l ( ζ , τ ) } − τ τ * 2 40 { [ 8 + 18 Z l 2 ( ζ , τ ) + 4 Z l 4 ( ζ , τ ) ] ierfc Z l ( ζ , τ ) − Z l ( ζ , τ ) [ 7 + 2 Z l 2 ( ζ , τ ) ] erfc Z l ( ζ , τ ) } − τ τ * 3 560 { 48 + 174 Z l 2 ( ζ , τ ) + 80 Z l 4 ( ζ , τ ) + 8 Z l 6 ( ζ , τ ) ierfc Z l ( ζ , τ ) − Z l ( ζ , τ ) [ 57 + 36 Z l 2 ( ζ , τ ) + 4 Z l 4 ( ζ , τ ) ] erfc Z l ( ζ , τ ) } 〉

(8)

T l , 4 * ( ζ , τ ) = 0 . 75 γ π τ τ * [ erfc Z l ( ζ , τ ) − 2 Z l ( ζ , τ ) ierfc Z l ( ζ , τ ) ]

(9)

T l , 5 * ( ζ , τ ) = 2 γ τ 〈 3 ierfc Z l ( ζ , τ ) − 2 τ * { 2 [ 1 + Z l 2 ( ζ , τ ) ] ierfc Z l ( ζ , τ ) − Z l ( ζ , τ ) erfc Z l ( ζ , τ ) } + 0 . 2 τ * 2 { [ 8 + 18 Z l 2 ( ζ , τ ) + 4 Z l 4 ( ζ , τ ) ] ierfc Z l ( ζ , τ ) − Z l ( ζ , τ ) [ 7 + 2 Z l 2 ( ζ , τ ) ] erfc Z l ( ζ , τ ) } 〉

(10)

T l , 6 * ( ζ , τ ) = 0 . 4 γ τ τ * 2 { [ 8 + 18 Z l 2 ( ζ , τ ) + 4 Z l 4 ( ζ , τ ) ] ierfc Z l ( ζ , τ ) − Z l ( ζ , τ ) [ 7 + 2 Z l 2 ( ζ , τ ) ] erfc Z l ( ζ , τ ) }

(11)

T l , 7 ∗ ( ζ , τ ) = 4 γ τ τ ∗ 〈 { 2 [ 1 + Z l 2 ( ζ , τ ) ] ierfc Z l ( ζ , τ ) − Z l ( ζ , τ ) erfc Z l ( ζ , τ ) } − 0 . 2 τ ∗ { [ 8 + 18 Z l 2 ( ζ , τ ) + 4 Z l 4 ( ζ , τ ) ] ierfc Z l ( ζ , τ ) − Z l ( ζ , τ ) [ 7 + 2 Z l 2 ( ζ , τ ) ] erfc Z l ( ζ , τ ) } 〉

(12)

T l , 8 ∗ ( ζ , τ ) = 0 . 8 γ τ 〈 3 ierfc Z l ( ζ , τ ) + τ ∗ { 2 [ 1 + Z l 2 ( ζ , τ ) ] ierfc Z l ( ζ , τ ) − Z l ( ζ , τ ) erfc Z l ( ζ , τ ) } − 0 . 4 τ ∗ 2 { [ 8 + 18 Z l 2 ( ζ , τ ) + 4 Z l 4 ( ζ , τ ) ] ierfc Z l ( ζ , τ ) − Z l ( ζ , τ ) [ 7 + 2 Z l 2 ( ζ , τ ) ] erfc Z l ( ζ , τ ) } 〉

(13)

T l , 9 ∗ ( ζ , τ ) = 0 . 8 γ τ τ ∗ 〈 3 { 2 [ 1 + Z l 2 ( ζ , τ ) ] ierfc Z l ( ζ , τ ) − Z l ( ζ , τ ) erfc Z l ( ζ , τ ) } − 0 . 4 τ ∗ { [ 8 + 18 Z l 2 ( ζ , τ ) + 4 Z l 4 ( ζ , τ ) ] ierfc Z l ( ζ , τ ) − Z l ( ζ , τ ) [ 7 + 2 Z l 2 ( ζ , τ ) ] erfc Z l ( ζ , τ ) } 〉

(14)

T l , 10 * ( ζ , τ ) = γ τ * 〈 3 π τ s [ erfc Z l ( ζ , τ ) − 2 Z l ( ζ , τ ) ierfc Z l ( ζ , τ ) ] − 4 τ { 2 [ 1 + Z l 2 ( ζ , τ ) ] ierfc Z l ( ζ , τ ) − Z l ( ζ , τ ) erfc Z l ( ζ , τ ) } 〉

(15)

Z 1 ( ζ , τ ) = ζ 2 τ , ζ ≥ 0 , Z 2 ( ζ , τ ) = − ζ 2 k * τ , ζ ≤ 0

(16)

ζ = z a 1 , τ = k 1 t a 1 2 , τ s = k 1 t s a 1 2 , τ * = τ τ s , γ = k * k * + K * , K * = K 2 K 1 , k * = k 2 k 1 , T 0 = q 0 a 1 K 1

(17)

a l = 3 k l t s , l = 1 , 2 are the effective depths of the friction couple elements heating; ⁴ ierfc ( x ) = e − x 2 / π − x erfc ( x ) , erfc ( x ) = 1 − erf ( x ) , where erf ( x ) is the Gauss error function. ²⁴

Thermal stresses

Spatio-temporal distributions of the quasi-static thermal stresses in the friction elements due to heterogeneous temperature fields T l , i ( z , t ) (equations (5)–(17)) will be found based on the theory of thermal bending of a thick plate with free edges.^25,26 In accordance with this theory, in every point of the strips 0 ≤ z ≤ a 1 —in a pad—and − a 2 ≤ z ≤ 0 —in a disk—in the directions Ox and Oy, the same normal thermal stresses act

σ l , i ( z , t ) = σ l , 0 σ l , i * ( ζ , τ ) , − a 2 ≤ z ≤ a 1 , 0 ≤ t ≤ t s

(18)

where

σ l , 0 = α l E l T 0 1 − ν l

(19)

σ l , i * ( ζ , τ ) = ε l , i * ( ζ , τ ) − T l , i * ( ζ , τ ) , a 2 * ≤ ζ ≤ a 1 * , 0 ≤ τ ≤ τ s

(20)

ε l , i * ( ζ , τ ) = [ ( − 1 ) l + 1 4 − 6 ζ ] M l , i ( 0 ) ( τ ) + [ ( − 1 ) l + 1 12 ζ − 6 ] M l , i ( 1 ) ( τ )

(21)

M l , i ( m ) ( τ ) = ∫ 0 a l * ζ m T l , i * ( ζ , τ ) d ζ , a l * = ( − 1 ) l + 1 a l / a 1 , i = 1 , 10 ¯ ; l = 1 , 2 ; m = 0 , 1

(22)

Substituting the dimensionless temperatures T l , i * ( ζ , τ ) (equations (6)–(17)) into the integrals in equation (22), we find

M l , 1 ( m ) ( τ ) = 4 3 γ τ τ * X l m + 1 ( τ ) { [ 3 / τ * − 2 ] I l ( m ) ( τ ) − 2 I l ( m + 2 ) ( τ ) + J l ( m + 1 ) ( τ ) }

(23)

M l , 2 ( m ) ( τ ) = 4 3 γ τ τ * X l m + 1 ( τ ) { 2 [ I l ( m ) ( τ ) + I l ( m + 2 ) ( τ ) ] − J l ( m + 1 ) ( τ ) }

(24)

M l , 3 ( m ) ( τ ) = γ X l m + 1 ( τ ) 〈 3 τ I l ( m ) ( τ ) − 0 . 5 τ τ * { 2 [ I l ( m ) ( τ ) + I l ( m + 2 ) ( τ ) ] − J l ( m + 1 ) ( τ ) } − τ τ * 2 40 { 8 I l ( m ) ( τ ) + 18 I l ( m + 2 ) ( τ ) + 4 I l ( m + 4 ) ( τ ) − 7 J l ( m + 1 ) ( τ ) − 2 J l ( m + 3 ) ( τ ) } − τ τ * 3 560 { 48 I l ( m ) ( τ ) + 174 I l ( m + 2 ) ( τ ) + 80 I l ( m + 4 ) ( τ ) + 8 I l ( m + 6 ) ( τ ) − 57 J l ( m + 1 ) ( τ ) − 36 J l ( m + 3 ) ( τ ) − 4 J l ( m + 5 ) ( τ ) } 〉

(25)

M l , 4 ( m ) ( τ ) = 0 . 75 π τ s τ * X l m + 1 ( τ ) [ J l ( m ) ( τ ) − 2 I l ( m + 1 ) ( τ ) ]

(26)

M l , 5 ( m ) ( τ ) = 2 γ τ X l m + 1 ( τ ) { 3 I l ( m ) ( τ ) − 2 τ * [ 2 I l ( m ) ( τ ) + 2 I l ( m + 2 ) ( τ ) − J l ( m + 1 ) ( τ ) ] + 0 . 2 τ * 2 [ 8 I l ( m ) ( τ ) + 18 I l ( m + 2 ) ( τ ) + 4 I l ( m + 4 ) ( τ ) − 7 J l ( m + 1 ) ( τ ) − 2 J l ( m + 3 ) ( τ ) ] }

(27)

M l , 6 ( m ) ( τ ) = 0 . 4 τ τ * 2 X l m + 1 ( τ ) [ 8 I l ( m ) ( τ ) + 18 I l ( m + 2 ) ( τ ) + 4 I l ( m + 4 ) ( τ ) − 7 J l ( m + 1 ) ( τ ) − 2 J l ( m + 3 ) ( τ ) ]

(28)

M l , 7 ( m ) ( τ ) = 4 γ τ τ * X l m + 1 ( τ ) { [ 2 I l ( m ) ( τ ) + 2 I l ( m + 2 ) ( τ ) − J l ( m + 1 ) ( τ ) ] − 0 . 2 τ * [ 8 I l ( m ) ( τ ) + 18 I l ( m + 2 ) ( τ ) + 4 I l ( m + 4 ) ( τ ) − 7 J l ( m + 1 ) ( τ ) − 2 J l ( m + 3 ) ( τ ) ] }

(29)

M l , 8 ( m ) ( τ ) = 0 . 8 τ X l m + 1 ( τ ) { 3 I l ( m ) ( τ ) + τ * [ 2 I l ( m ) ( τ ) + 2 I l ( m + 2 ) ( τ ) − J l ( m + 1 ) ( τ ) ] − 0 . 4 τ * 2 [ 8 I l ( m ) ( τ ) + 18 I l ( m + 2 ) ( τ ) + 4 I l ( m + 4 ) ( τ ) − 7 J l ( m + 1 ) ( τ ) − 2 J l ( m + 3 ) ( τ ) ] }

(30)

M l , 9 ( m ) ( τ ) = 0 . 8 τ τ * X l m + 1 ( τ ) { 3 [ 2 I l ( m ) ( τ ) + 2 I l ( m + 2 ) ( τ ) − J l ( m + 1 ) ( τ ) ] − 0 . 4 τ * [ 8 I l ( m ) ( τ ) + 18 I l ( m + 2 ) ( τ ) + 4 I l ( m + 4 ) ( τ ) − 7 J l ( m + 1 ) ( τ ) − 2 J l ( m + 3 ) ( τ ) ] }

(31)

M l , 10 ( m ) ( τ ) = γ τ * X l m + 1 ( τ ) { 3 π τ s [ J l ( m ) ( τ ) − 2 I l ( m + 1 ) ( τ ) ] − 4 τ [ 2 I l ( m ) ( τ ) + 2 I l ( m + 2 ) ( τ ) − J l ( m + 1 ) ( τ ) ] }

(32)

where

I l ( k ) ( τ ) = ∫ 0 Y l ( τ ) x k ierfc ( x ) dx = π − 1 / 2 L l ( k ) ( τ ) − J l ( k + 1 ) ( τ ) , k = 0 , 7 ¯

(33)

J l ( k ) ( τ ) = ∫ 0 Y l ( τ ) x k erfc ( x ) dx , k = 0 , 8 ¯

(34)

L l ( k ) ( τ ) = ∫ 0 Y l ( τ ) x k e − x 2 dx , k = 0 , 7 ¯

(35)

X l ( τ ) = 1 Z l ( 1 , τ ) , Y l ( τ ) = a l * X l ( τ )

(36)

Function Z l ( 1 , τ ) , l = 1 , 2 is determined from relation (16). Using the recursive formulas^27,28

J l ( k ) ( τ ) = 1 ( k + 1 ) { k ( k − 1 ) 2 J l ( k − 2 ) ( τ ) + [ Y l 2 ( τ ) − k 2 ] Y l k − 1 ( τ ) erfc Y l ( τ ) − 1 π Y l k ( τ ) e − Y l 2 ( τ ) }

(37)

L l ( k ) ( τ ) = 0 . 5 [ ( k − 1 ) L l ( k − 2 ) ( τ ) − Y l k − 1 ( τ ) e − Y l 2 ( τ ) ] , k = 2 , 3 , …

(38)

J l ( 0 ) ( τ ) = π − 1 / 2 − ierfc Y l ( τ ) , J l ( 1 ) ( τ ) = 0 . 25 erf Y l ( τ ) − 0 . 5 Y l ( τ ) ierfc Y l ( τ )

(39)

L l ( 0 ) ( τ ) = 0 . 5 π [ 1 − erfc Y l ( τ ) ] , L l ( 1 ) ( τ ) = 0 . 5 [ 1 − e − Y l 2 ( τ ) ]

(40)

The integrals J l ( k ) ( τ ) (equation (34)) and L l ( k ) ( τ ) (equation (35)) are written in the form

J l ( 2 ) ( τ ) = 3 − 1 〈 Y l 3 ( τ ) erfc Y l ( τ ) + π − 1 / 2 { 1 − [ Y l 2 ( τ ) + 1 ] e − Y l 2 ( τ ) } 〉

(41)

J l ( 3 ) ( τ ) = 4 − 1 { Y l 4 ( τ ) erfc Y l ( τ ) + 0 . 75 erf Y l ( τ ) − π − 1 / 2 Y l ( τ ) [ Y l 2 ( τ ) + 1 . 5 ] e − Y l 2 ( τ ) }

(42)

J l ( 4 ) ( τ ) = 5 − 1 〈 Y l 5 ( τ ) erfc Y l ( τ ) + π − 1 / 2 { 2 − [ Y l 4 ( τ ) + 2 Y l 2 ( τ ) + 2 ] e − Y l 2 ( τ ) } 〉

(43)

J l ( 5 ) ( τ ) = 6 − 1 { Y l 6 ( τ ) erfc Y l ( τ ) + 1 . 875 erf Y l ( τ ) − π − 1 / 2 Y l ( τ ) [ Y l 4 ( τ ) + 2 . 5 Y l 2 ( τ ) + 3 . 75 ] e − Y l 2 ( τ ) }

(44)

J l ( 6 ) ( τ ) = 7 − 1 〈 Y l 7 ( τ ) erfc Y l ( τ ) + π − 1 / 2 { 6 − [ Y l 6 ( τ ) + 3 Y l 4 ( τ ) + 6 Y l 2 ( τ ) + 6 ] e − Y l 2 ( τ ) } 〉

(45)

J l ( 7 ) ( τ ) = 8 − 1 { Y l 8 ( τ ) erfc Y l ( τ ) + 6 . 5625 erf Y l ( τ ) − π − 1 / 2 Y l ( τ ) [ Y l 6 ( τ ) + 3 . 5 Y l 4 ( τ ) + 8 . 75 Y l 2 ( τ ) + 13 . 125 ] e − Y l 2 ( τ ) }

(46)

J l ( 8 ) ( τ ) = 9 − 1 〈 Y l 9 ( τ ) erfc Y ( τ ) l + π − 1 / 2 { 24 − [ Y l 8 ( τ ) + 4 Y l 6 ( τ ) + 12 Y l 4 ( τ ) + 24 Y l 2 ( τ ) + 24 ] e − Y l 2 ( τ ) } 〉

(47)

L l ( 2 ) ( τ ) = 0 . 25 π erf Y l ( τ ) − 0 . 5 Y l ( τ ) e − Y l 2 ( τ ) , L l ( 3 ) ( τ ) = 0 . 5 { 1 − [ 1 + Y l 2 ( τ ) ] e − Y l 2 ( τ ) }

(48)

L l ( 4 ) ( τ ) = 0 . 375 π erf Y l ( τ ) − 0 . 5 Y l ( τ ) [ Y l 2 ( τ ) + 1 . 5 ] e − Y l 2 ( τ )

(49)

L l ( 5 ) ( τ ) = 1 − [ 0 . 5 Y l 4 ( τ ) + Y l 2 ( τ ) + 1 ] e − Y l 2 ( τ )

(50)

L l ( 6 ) ( τ ) = 0 . 9375 π erf Y l ( τ ) − Y l ( τ ) [ 0 . 5 Y l 4 ( τ ) + 1 . 25 Y l 2 ( τ ) + 1 . 875 ] e − Y l 2 ( τ )

(51)

L l ( 7 ) ( τ ) = 3 − [ 0 . 5 Y l 6 ( τ ) + 1 . 5 Y l 4 ( τ ) + 3 Y l 2 ( τ ) + 3 ] e − Y l 2 ( τ )

(52)

Taking into account the functions J l ( k ) ( τ ) and L l ( k ) ( τ ) (equations (39)–(52)), from relation (33), the integrals were determined I l ( k ) ( τ ) , l = 1 , 2 ; k = 0 , 7 ¯

I l ( 0 ) ( τ ) = 0 . 5 π − 1 / 2 Y l ( τ ) e − Y l 2 ( τ ) + 0 . 25 erf Y l ( τ ) − 0 . 5 Y l 2 ( τ ) erfc Y l ( τ )

(53)

I l ( 1 ) ( τ ) = 6 − 1 〈 π − 1 / 2 { 1 + [ 2 Y l 2 ( τ ) − 1 ] e − Y l 2 ( τ ) } − 2 Y l 3 ( τ ) erfc Y l ( τ ) 〉

(54)

I l ( 2 ) ( τ ) = 8 − 1 { π − 1 / 2 Y l ( τ ) [ 2 Y l 2 ( τ ) − 1 ] e − Y l 2 ( τ ) + 0 . 5 erf Y l ( τ ) − 2 Y l 4 ( τ ) erfc Y l ( τ ) }

(55)

I l ( 3 ) ( τ ) = 5 − 1 〈 0 . 5 π − 1 / 2 { 1 + [ 2 Y l 4 ( τ ) − Y l 2 ( τ ) − 1 ] e − Y l 2 ( τ ) } − Y l 5 ( τ ) erfc Y l ( τ ) 〉

(56)