- Email: [email protected]
Laser pulse heating of steel surface and flexural wave analysis
Optics and Lasers in Engineering 37 (2002) 63–83
Laser pulse heating of steel surface and flexural wave analysis B.S. Yilbas*, M. Faisal, S.Z. Shuja, A.F.M. Arif King Fahd University of Petroleum and Minerals, Dhahran, Saudi Arabia Received 13 July 2000; received in revised form 1 June 2001; accepted 1 August 2001
Abstract Laser gas-assisted material processing finds wide application in industry. The modelling of heating, elastic response of the substrate material, and the wave analysis gives insight into the laser workpiece interaction. In the present study, laser gas-assisted heating of steel is considered. The normal component of the thermal stress is taken as the source of load for the flexural wave generation in the material. The flexural wave generated is simulated and the wave characteristics are analyzed at four locations at the workpiece surface. The numerical scheme employing a control volume approach is introduced when solving the governing equations of flow and heat transfer while finite element and spectran element methods are used when solving the stress and wave equations. It is found that the normal component of the stress is tensile. The dispersion effect of the workpiece material, interference of the reflected beam, and partial overlapping of second mode of the travelling wave enable to identify a unique pattern in the travelling wave in the substrate. r 2001 Elsevier Science Ltd. All rights reserved. Keyword: Laser heating flexural wave generation
1. Introduction Lasers are widely used in industry due to their precision of operation, high speed of processing, and low cost. When the laser beam is absorbed by the substrate material, temperature of the irradiated region rises, which in turn results in thermal stress generation due to thermal expansion of the workpiece material. Thermal loading produces flexural waves propagating in the substrate. The workpiece acts as *Corresponding author. Tel.: +966-3-860-2540; fax: +966-3-860-2949. E-mail address: [email protected] (B.S. Yilbas). 0143-8166/02/$ - see front matter r 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 1 4 3 - 8 1 6 6 ( 0 1 ) 0 0 1 4 9 - X
64
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
Nomenclature a Gaussian parameter A cross-sectional area of workpiece, m2 cp specific heat capacity, J=kg K C1 ; C2 ; Cm coefficients in the k–e turbulence model d diameter of the nozzle, m D elasticity matrix Dp surface displacement, m f1 ; f2 ; fm coefficients in the low Reynolds number, k–e model Df frequency step, Hz E modulus of elasticity, MPa G rate of generation of k; W=m3 I0 peak power intensity, W=m2 I moment of inertia of cross section, m4 k turbulent kinetic energy, W=m3 k wave number (Eqs. 26–28) M bending moment, Nm N number of data points in power spectrum p pressure, Pa Pr varible Prandtl number (function of temperature) q heat flux, W=m2 Re Reynolds number r distance in the radial direction, m S unsteady spatially varying source (Eq. (9)), W=m3 t time, s T temperature, K U arbitrary velocity, m=s V radial velocity, m=s W axial velocity, m=s x arbitrary direction, m zn distance to the nearest wall, m z distance in the axial direction, m DTðr; z; tÞ Temperature rise at a point ðr; zÞ at time¼ t with respect to that at t ¼ 0 Tref Reference temperature at t ¼ 0 ff B g the applied body force fPg the applied pressure vector % fFg concentrated nodal forces to the element fdUg virtual displacement fdUs g virtual displacement on the boundary where pressure is prescribed % fdUg virtual displacement of boundary nodes where concentrated load is prescribed ½B strain displacement gradient matrix % fUg nodal displacement vector
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
65
½N matrix of shape (or interpolation) functions R fF th g element thermal load vector, 8 ½BTR½D½eth d8 T B fFgb element applied body force 8 ½N ½f d8 R vector, s T fFg element pressure vector, R ½Nn ½P d ½Ke element stiffness matrix, 8 ½BT ½D½B d8 Greek letters s feg feth g e0 a u G d e Z l k k m
r F f f o
stress vector, MPa total strain vector thermal strain vector equivalent strain coefficient of thermal expansion, 1/K Poisson’s ratio diffusion coefficient, m2 =s absorption coefficient, 1=m energy dissipation, W=kg damping coefficient thermal conductivity, W=mK turbulence intensity shear correction factor (Eqs. 22; 23) variable dynamic viscosity (function of temperature for laminar), N s=m2 effective viscosity, ðml þ mt Þ; N s=m2 variable kinematic viscosity (function of temperature for laminar), m2 =s density (function of temperature and pressure for gas), kg=m3 viscous dissipation, W=m3 arbitrary variable bending slope circular frequency, Hz
Subscript amb i; j jet l max t
ambient arbitrary direction conditions at jet inlet laminar maximum turbulent
me n
a dispersive media influencing the wave propagation. The wave response of the workpiece is governed by three system parameters namely, damping, mass, and stiffness. Consequently, the geometry, the mechanical properties, and the position of the workpiece determine the wave propagation characteristics.
66
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
In order to analyze the laser-induced flexural waves, the stress field generated in the laser-irradiated region needs to be analyzed first. This is because the load generating the flexural wave is due to the axial and transverse stress components in the axisymmetric system without torsional motion. Considerable research studies were conducted on the laser-induced thermal stresses in solids. The laser-induced thermal stresses on steel surfaces were investigated by Yilbas et al. [1]. They showed that excessive stress levels occurred in the irradiated spot, which in turn resulted in the deformation at the surface of the substrate material. The thermal stresses generated during the laser quenching was investigated by Wang et al. [2]. They indicated that the residual stress developed in the surface vicinity resulted in the hardening of the surface. The experimental and numerical studies of fracture initiation in thin aluminum oxide ceramics during laser cutting were carried out by Li and Sheng [3]. They introduced plane stress model when predicting the stresses along the cut edges. They showed that low cutting speed reduces the stress level and minimizes the possible fracture initiation at the cut edges. An analytical solution of a two-dimensional thin coating attached to a substrate was obtained by Elperin and Rudin [4]. They indicated that the temperature difference across the coating generated excessive stress levels, which in turn caused inelastic deformation of the coating. Modest [5] investigated the elastic and viscoelastic stresses during laser drilling of ceramics. He indicated that the material softened in the region closer to the melt zone. The temperature field generated in the laser-irradiated region plays an important role in stress development and flexural wave initiation in the substrate. Considerable research studies were carried out to investigate the temperature field in the substrate during laser heating process. Heat conduction losses in laser cutting of metals was investigated by Schultz et al. [6]. They showed that the theoretical predictions agreed well with the experimental findings. Shuja and Yilbas [7] investigated the gas-assisted laser heating of steel substrate. They introduced the three-dimensional conjugate heating model and predicted the temperature field in both the solid and the gas side. They indicated that the impinging gas jet had little influence on the surface temperature rise. Neto Diniz and Lima [8] employed a three-dimensional parabolic heat diffusion equation to predict the temperature field during laser pulse heating of steel surfaces. They indicated that the oscillation in the pulse had no significant effect on the resulting heat-affected zone depth. Yilbas [9] introduced a three-dimensional heating model with a stationary Gaussian laser beam source. The phase change process was employed in the analysis. He indicated that the temperature below the surface became higher than that of the surface as the rate of surface evaporation increased. The laser-induced waves in solids received considerable attention, since it served as a nondestructive tool for materials evaluation [10,11]. Several analytical methods were introduced when modelling the thermoelastic process that took place during the laser heating pulse. Some of these include Green function method, and Laplace and Henkel transformations [12,13]. However, the thermal stresses developed in the substrate vary with time and it is, in general, difficult to present in a simple mathematical function. Therefore, the analytical solution to the transient problem
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
67
with nonordered input load becomes very difficult to obtain. Consequently, the general trend is to solve the governing equations numerically. Dubois et al. [14] investigated numerically the laser thermoelastic generation of ultrasound in an orthotropic medium. They compared the predictions with the experimental results and indicated that their predictions agreed well with the experimental findings. They further argued that the model developed could be used to optimize the laser generation of ultrasound in various materials. Cheng et al. [15] studied the dynamic thermoelastic response of pulsed photothermal deformation deflection detection for Q-switched laser pulses. They indicated that when the laser pulse rise time decreases in the order of 10 ns or less, the pulsed photothermal deformation deflection signal reflected a totally dynamic wave behavior. The laser-induced waves such as Lamb waves can also be used to determine the properties of the substrate material. The Lamb wave is generated due to normal and shear loads applied to the surface of the half-space [16]. Dewhurst et al. [17] introduced a Lamb wave technique to describe the measurement of the film thickness on the substrate material. A theoretical study for modelling the thermoelastic excitation of transient Lamb waves propagating along the principal direction in an orthotropic plate was carried out by Cheng and Berthelod [18]. They provided a quantitative analysis for noncontact and nondestructive detection of the elastic stiffness properties of the machine-made paper by using the laser-generated Lamb wave technique. In the light of the above arguments, the present study is carried out to investigate the flexural wave propagation in the substrate during the laser heating process. The thermal stresses developed in the substrate are predicted first and the wave generated due to axial stress component in the surface vicinity of the substrate is analyzed using a spectral finite element method. Since the workpiece has a certain length, the reflecting waves from the free supported ends of the workpiece are also considered in the analysis.
2. The mathematical model The gas-assisted laser pulse heating process is considered when modelling the laser workpiece interaction provided that conduction limited heating case is taken into account in the present study [19]. Fig. 1 shows the schematic view of the heating process. 2.1. Flow and heat conduction equations The unsteady flow equations need to be solved to obtain the flow field due to axisymmetric gas jet impingement. The continuity and momentum equations are qr q þ ðrUi Þ ¼ 0 qt qxi
ð1Þ
68
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
Nozzle
r Nozzle exit 1 mm dia. Air gas jet velocity100 m/s Stand off 2 mm
z
Solid gas interface
Solid thickness 2 mm Stress free boundary Free support
Stress free boundary Free support
Fig. 1. A view of impinging gas and laser workpiece arrangement.
and
q q qp q qUj ðrUj Þ þ ðrUi Uj Þ ¼ þ ðm þ mÞ ; qt qxi qxj qxi t qxi
ð2Þ
where mt is the eddy viscosity which has to be specified by a turbulence model. In the present study, low Reynolds number k–e turbulence model is introduced to account for the turbulence effect of the impinging gas, which is air [19]. The partial differential equation governing the transport of thermal energy has the form qT q q mt m qT þ ðUi TÞ ¼ þ : ð3Þ qt qxi qxi Prt Pr qxi
2.1.1. Boundary conditions for flow equations Laminar boundary conditions are set for the mean-flow variables, and the boundary conditions k ¼ 0 and de=dz ¼ 0 are applied at the wall. The low-Reynolds number extension does not employ wall functions; therefore, the grid employed normal to the main flow direction needs to be distributed so as to give a high concentration of grid cells near the wall, with the wall-adjacent node positioned at zþ ¼ rzun =mp1:0: Inlet to control volume. Ui ¼ specified
and
T ¼ constant:
ð4Þ
The kinetic energy of turbulence is estimated according to some fraction of the square of the average inlet velocity [20] k ¼ lu%2 ; where u% is the average inlet velocity and l is a fraction.
ð5Þ
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
69
The dissipation is calculated according to the equation [11] k3=2 ; ð6Þ bd where d is the inlet diameter. The values l ¼ 0:03 and b ¼ 0:005 are commonly used and may vary slightly in the literature [20]. Outlet to control volume It is assumed that the flow extends over a sufficiently long domain; therefore, it is fully developed at the exit section. Thus, for any variable f the condition is qðrfÞ ¼ 0; ð7Þ qx where x is the arbitrary outlet direction. Symmetry axis The radial derivative of the variables is set to zero at the symmetry axis, i.e. qf ¼ 0 and V ¼ 0: ð8Þ qr Solid fluid interface The temperature at the solid–gas interface is considered as the same, i.e. qTwgas qTwsolid Twsolid ¼ Twgas and Ksolid ¼ Kgas : qz qz The unsteady heat conduction equation. q q qT ðcp rTÞ ¼ K þ S; ð9Þ qt qxi qxi e ¼ Cm
where S is the unsteady spatially varying laser output power intensity distribution and is considered as Gaussian with 1=e2 points equal to 0:375 mm; from the center of the beam. Therefore, S is stated as follows: 2 I0 r S ¼ pffiffiffiffiffiffiffiffi exp 2 d expð dzÞf ðtÞ; ð10Þ a 2pa pffiffiffiffiffiffiffiffi where ðI0 = 2paÞ expðr2 =a2 Þ is the intensity distribution across the surface, expð dzÞ is the absorption function, and f ðtÞ is the function accommodating the time variation of the pulse shape. The pulse properties are given in Table 1. Boundary conditions for the heat conduction equation Convection with a constant coefficient for still air is considered at the z ¼ zmax boundary for the plate. The continuity of temperature between the solid and the gas is enforced at the interface;
Table 1 Pulse properties used in the simulation Rise time (ms)
Fall time (ms)
Total pulse length (ms)
Peak power intensity ðW=m2 Þ
0.2
1.08
1.48
109
70
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
and a constant temperature, T ¼ Tamb ; is assumed for the distance far away from the laser source. 2.1.2. Variable properties An equation of state is used for the impinging gas and the specific heat capacity and thermal conductivity for both air and steel were considered only as a function of temperature. The temperature dependence of properties is given in [21]. 2.1.3. Calculation of the flow field variables The control volume approach is used when solving the governing equations of flow and heat transfer numerically. The differential equation is integrated over the control volume to yield the discretization equation. The main reasons for choosing the control-volume formulation for flow field are its simplicity and easy physical interpretation [22]. The discretization process is not given here due to lengthy arguments, but refer to [23]. The grid used in the present calculations has 38 70 mesh points, provided that the grid independent test is satisfied. The details of the grid orientation and grid independent test results are given in [9]. The two problems of determining the pressure and satisfying continuity are overcome by adjusting the pressure field so as to satisfy continuity. This arrangement gives a convenient way of handling the pressure linkages through the continuity equation and is known as the SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algorithm. The details of this algorithm are found in [24]. The governing equation for heat conduction in solid (Eq. (9)) can be written in the form of flow equations. Thus, the discretization procedure leads to algebraic equations of the form similar to flow equations with temperature T replacing the general variable. 2.2. Thermal stress modelling The temperature field results in the thermal stresses being generated in the substrate, which can lead to the elastic–plastic displacement in the substrate material. The stress is related to strains by fsg ¼ ½Dfee g;
ð11Þ
where fsg is the stress vector, and ½D is the elasticity matrix. fee g ¼ feg feth g; where feg is the total strain vector and feth g is the thermal strain vector. Eq. (11) may also be written as feg ¼ ½D 1 fsg þ feth g:
ð12Þ
Since the present case is axially symmetric, and the material is assumed to be isotropic, the above stress–strain relations can be written in cylindrical
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
coordinates as 1 err ¼ ½srr uðsyy þ szz Þ þ aDTðr; z; tÞ; E 1 eyy ¼ ½syy uðsrr þ szz Þ þ aDTðr; z; tÞ; E 1 ezz ¼ ½szz uðsrr þ syy Þ þ aDTðr; z; tÞ; E 1 erz ¼ srz : G
71
ð13Þ
where E; u; and a are the modulus of elasticity, Poisson’s ratio, and coefficient of thermal expansion, respectively. DTðr; z; tÞ represents the temperature rise at a point ðr; zÞ at time ¼ t with respect to that at t ¼ 0 corresponding to a stress-free condition. A typical component of thermal strain from Eq. (13) is eth ¼ aDTðr; z; tÞ ¼ aðTðr; z; tÞ Tref Þ; where Tref is the reference temperature at t ¼ 0: When a is a function of temperature, then Eq. (14) becomes Z T th e ¼ aðTÞ dT:
ð14Þ
ð15Þ
Tref
The present study uses a mean or weighted-average value of a; such that eth ¼ a% ðTÞðT Tref Þ;
ð16Þ
where a% ðTÞ is the mean value of coefficient of thermal expansion and is given by RT Tref aðTÞ dT a% ðTÞ ¼ : ð17Þ Tðr; z; tÞ Tref To develop a finite element procedure for stress computation, the standard displacement-based finite element method is used. The basis of this approach is the principle of virtual work, which states that the equilibrium of any body under loading requires that for any compatible small virtual displacements (which are zero at the boundary points and surfaces and correspond to the components of displacements that are prescribed at those points and surfaces) imposed on the body in its state of equilibrium, the total internal virtual work or strain energy ðdUÞ is equal to the total external work due to the applied thermally induced loads ðdVÞ; i.e. dU ¼ dV: For the static analysis of problems having linear geometry and thermoelastic material behavior, one can derive the following equation using standard procedure [25]. Z ðfdegT ½Dfeg fdegT ½Dfeth gÞ d8 8 Z Z X % T fFg; % ¼ fdUgT ff B g d8 þ fdUs gT fPg dO þ fdUg ð18Þ 8
72
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
% the where ff B g is the applied body force, fPg the applied pressure vector, fFg concentrated nodal forces to the element, fdUg the virtual displacement, fdUs g the % the virtual displacement on the boundary where pressure is prescribed, and fdUg virtual displacement of boundary nodes where concentrated load is prescribed. The strains may be related to the nodal displacement by % feg ¼ ½BfUg;
ð19Þ
% the nodal displacement where ½B is the strain displacement gradient matrix, and fUg vector. The displacements within the elements are related to the nodal displacement by % fUg ¼ ½NfUg; ð20Þ where ½N is the matrix of shape (or interpolation) functions. Eq. (18) can be reduced to the following matrix form: % fF th g ¼ fFgb þ fFgs þ fFg % ½Ke fUg ð21Þ R R T T where ½Ke ¼ 8 ½B ½D½B d8 is the element stiffness fF th g ¼ 8 ½B ½D R matrix, b T B th ½A d8 the element Rthermal load vector, fFg ¼ 8 ½N ½f d8 the element body force vector, fFgs ¼ ½Nn T ½P d the element pressure vector, and ½Nn ¼ matrix of shape functions for normal displacement at the boundary surface. Assembly of element matrices and vectors of Eq. (21) yields % ¼ fRg; % ½Kfdg % and fRg % are the global stiffness matrix, global nodal displacement where ½K; fdg vector, and global nodal load vector, respectively. Solution of the above set of simultaneous algebraic equations gives unknown nodal displacements and reaction forces. 2.3. Wave analysis In the analysis, the workpiece is divided into a number of elements and each element can be treated as uniform. The displacement vector is then transformed from time domain into frequency domain using fast Fourier transformation (FFT). The transform response is fed back into the governing differential equations as solved at each frequency component using finite element method (FEM) and the transient response is recovered by using the inverse FFT. 2.3.1. Formulation The solution in the wave problems is general function of space and time. If the time variation of the solution is focussed on a particular point in space, then it can be presented by a spectral function. The spectral formulation starts with the equations of motion of the workpiece (as free supported beam) including inertia terms. After dividing the workpiece into a number of elements, cross sectional area and second moment of area in the center of each element can be considered as constant
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
73
over the entire element. Therefore, the equation of motion can be written as [26] 2 q y qf q2 y GAk ð22Þ ¼ rA 2 q 2 qx qx qt and q2 f qy q2 f f ¼ rI þ GAk : qx2 qx qt2
ð23Þ
The displacement and the bending slope can be written in the spectral form as [26] X ; yðx; tÞ ¼ yn ðx; on Þeion t ð24Þ n
and fðx; tÞ ¼
X
;
fn ðx; on Þeion t :
ð25Þ
n
The spectral components yn ; and fn ; have the solutions ;
yn ¼ A1 e ik1 x þ B1 e ik2 x þ C1 eik1 x þ D1 eik2 x
ð26Þ
and ;
fn ¼ A2 e ik1 x þ B2 e ik2 x þ C2 eik1 x þ D2 eik2 x ;
ð27Þ
where k1 and k2 are the wave numbers and A1 ; A2 ; B1 ; B2 ; C1 ; C2 and D1 ; D2 are the frequency-dependent coefficients, which are complex in nature. The first term in Eqs. (26) and (27) represents the waves moving in the forward direction while the last two describe the backward moving waves. Substituting yn ; and fn ; into Eqs. (22) and (23), and considering the load to be applied at the modes, yields GAkð k21 A1 þ ik1 A2 Þ ¼ o2 rAð A1 Þ; GAkð k21 B1 þ ik1 B2 Þ ¼ o2 rAð B1 Þ; GAkð k21 C1 þ ik1 C2 Þ ¼ o2 rAð C1 Þ; GAkð k21 D1 þ ik1 D2 Þ ¼ o2 rAð D1 Þ; where the coefficients are as follows:
A2 ¼
o2 rA þ k21 GAk A1 ; ik1 GAk
B2 ¼
o2 rA þ k21 GAk B1 ; ik1 GAk
ð28Þ
74
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
2 o rA k21 GAk C1 ; ik1 GAk 2 o rA k21 GAk D2 ¼ D1 : ik1 GAk
C2 ¼
The shear force and momentum equations can be written as 2
;
d f Q’ ¼ EI 2 o2 rI f dx
ð29Þ
and ;
’ ¼ EI d f : M dx
ð30Þ
The dynamic stiffness is obtained by first relating the coefficients to the nodal displacement as ;
; ;
;
;
½A; B; C; D ¼ a ½y1 ; f1 ; y2 ; f2 ;
ð31Þ
where ½a; ¼ ½b and ½b 1 are given in [27]. The nodal loads are obtained by using the shear force and momentum equations Eqs. (29) and (30), i.e. ;
;
½F T ¼ ½k ½dT ; where ;
;
;
;
;
½F T ¼ ½Q1 M1 y Qn Mn and ;
;
;
;
½dT ¼ ½y1 f1 y y1 f1 where n is the number of nodes. 2.3.2. Initial and boundary conditions Initially the displacements and velocities are set to zero. The free supported boundary condition at the workpiece ends is considered. In this case, at free ends of the workpiece ðx ¼ 0 ¼ 1Þ; the bending moment is set to zero at all frequencies, i.e. d f; ðI; oÞ ¼ 0: dx In order to observe the effect of reflecting waves from the free ends of the workpiece, one end of the workpiece is assumed to be at infinity. This boundary condition (element extends to infinity) requires that the element behaves as a throwoff beam, which acts as a conduit for energy out of the system.
75
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
2.3.3. Method of solution In order to obtain the transient response of the workpiece using a spectral finite element method, the input excitation signal, which is the load generated due to normal component of the thermal stress, at discrete time steps needs to be converted into frequency domain by FFT. The sampling rate Dt must be in accordance with Nyquist frequency of the signal. The frequency step is obtained from 1 Df ¼ ; NDt where N is the number of data points in the power spectrum. The signal is bondlimited between the minimum frequency ðfmin Þ and the maximum frequency ðfmax Þ: Beginning from the minimum frequency, the assembled dynamic stiffness matrix and the assembled load vector in frequency domain are obtained by the finite element method, i.e.: ;
; ;
;
;
½F ¼ ½K ½d
ð32Þ
or ;
½d ¼ ½K 1 ½F : After applying the boundary conditions to Eq. (32), ½K is inverted by Gauss reduction method. Consequently, the response is obtained at fmin : This procedure is repeated for other frequencies upto fmax : Hence, the displacement vector ½d ; is obtained at the desired element. The frequency-dependent response is reconstructed with time domain by using inverse Fourier transformation method. In order to improve the stability of the computation, damping is introduced through the wave number k in the form k ¼ k0 ð1 iZÞ where k0 is the undamped value and Z is taken as 0.025. The material properties are given in Table 2 while workpiece size and number of elements used in the simulation are given in Table 3. Table 2 Properties of steel used in the simulation and the number of elements employed E (GPa)
G (GPa)
r ðkg=m3 Þ
k
u
Z
207
77.6
7836
0.67
0.3
0.025
Table 3 Workpiece size and number of elements used in simulations Length (m)
Width (m)
Thickness (m)
Number of elements
1
0.2
0.002
8
76
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
3. Results and discussions The thermal stresses inside the substrate due to the temperature field developed by a laser heating pulse is modelled. The flexural wave generation in the substrate material due to the normal component of the thermal stress is analyzed. In order to investigate the wave propagation along the transverse direction in the workpiece, six locations with three located at each side from the center of the workpiece along the transverse direction are considered. This is necessary, since the workpiece material disperses the travelling waves in amplitude and frequency domain. The wave reflections from both ends of the workpiece also modify the travelling wave properties. Consequently, to analyze the influence of reflected waves emanating from the ends of the workpiece, a case in which one end of the workpiece is located at infinity is considered. This provides no reflection of the wave from one end of the workpiece. Fig. 2 shows the temporal variation of the surface temperature at the center of the heated spot (at r ¼ 0) for convective and nonconvective boundary conditions at the surface. The temperature rises rapidly in the early heating period and the rate of rise of surface temperature reduces as the heating progresses. In the cooling cycle, the decay rate of the surface temperature is faster in the beginning of the cooling cycle and the decay rate slows as the cooling period progresses. Consequently, the material response to the heating pulse changes in the cooling cycle as compared to that in the heating cycle. This is because once the laser source disappears, the rate of energy diffused from the surface vicinity of the substrate to the bulk through the conduction becomes large immediately after the laser pulse ends, i.e. the large temperature gradient attainment in the surface region during the heating cycle accelerates the conduction cooling of the surface vicinity. In the case of the long cooling period, the
Fig. 2. Temporal variation of surface temperature at the center of the irradiated spot with and without convective boundary condition at the surface.
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
77
temperature gradient in the surface vicinity is low and the energy diffusion from the surface vicinity to the bulk of the substrate slows down. When comparing temperature profiles with and without convective boundary condition at the surface, it is evident that both temperature profiles are identical. This indicates that convective cooling of the surface is negligibly smaller than the energy input to the substrate material due to the absorption of a laser beam as consistent with the previous work [9]. Fig. 3a shows the temporal variation of the normal (axial) component of the thermal stress while Fig. 3b shows the temporal variation of the equivalent strain in the surface vicinity. The equivalent strain increases as the heating progresses
25
STRESS (MPa)
18
11
4
-3
-10 0
0.5
1
1.5
2
1.2
1.6
TIME (ms)
(a) 0.025
STRAIN
0.020
0.015
0.010
0.005
0.000 0
(b)
0.4
0.8
TIME (ms)
Fig. 3. (a) Temporal variation of load generated in the surface vicinity of the workpiece due to normal component of the thermal stress; (b) Temporal variation of strain in the surface vicinity of the workpiece.
78
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
provided that the increase in strain is continuous even in the cooling cycle. It reduces rapidly as the cooling period progresses further. This is because the elastic response of the material is not as fast as the temperature field, i.e. the expansion is continuous for a while after the laser source vanishes. The normal component of the stress varies considerably with time. It results in two maxima; one in the heating cycle and the other in the cooling cycle. The normal component of the stress in the surface vicinity is tensile and it reaches as high as 20 MPa; which is less than the yield stress of the workpiece material. Fig. 4a shows the load variation with time while Fig. 4b shows the power spectrum density curve. The temporal variation of the load is similar to the normal stress component as shown in Fig. 3a, i.e. the value of the stress component is multiplied by the size of the irradiated spot to obtain the load. Since the total power in a signal is constant, the low and high frequency limits of the signal are determined from Fig. 4b. In the present case, Nyquist frequency (fc ¼ 1=2Dt; where Dt is the sampling rate) is selected as 2000. It should be noted that fc is the maximum frequency contained in the signal. Fig. 5 shows the flexural wave amplitude variation with time at four equally spaced locations at the workpiece surface. The first location is 12:5 cm away from the center of the heated spot and following locations are 12:5 cm away from the first location, respectively. The wave amplitude varies with time and no dominant pattern is observed in the wave. This is because the initial wave form at the center of the workpiece is modified as it travels along the workpiece due to: (i) the dispersion effect of the workpiece material, (ii) overlapping of wave modes, and (iii) the reflected waves from the free ends of the workpiece that interfere with the travelling wave. The group speed of the first wave mode is about 3125 m=s; therefore, the reflection from the free ends takes about 0:16 10 3 s to reach the center of the workpiece. Consequently, the reflected wave modifies the travelling wave in the time and amplitude domain. As the location from the first point changes to the following
Fig. 4. (a) Temporal variation of load generated due to normal component of the thermal stress; (b) Power spectrum density of thermal load.
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
79
Fig. 5. The wave generated at different locations of the workpiece.
Location number
1
2
3
4
Distance away from the irradiated spot center (m)
0.125
0.25
0.375
0.5
point along the surface, the amplitude of the wave decreases and the travelling wave pattern changes. This indicates that the travelling wave vanishes as it travels towards the free ends of the workpiece, i.e. in the present study, damping coefficient of the workpiece material is taken as 0.025. Moreover, at locations close to the free ends of the workpiece, the wave appears as in high frequency mode. This is because the travelling wave has a low amplitude in these locations and reflected waves modify easily the travelling wave pattern. In order to identify the influence of the reflected waves on the travelling wave, throw-off element is introduced at one end of the workpiece. This corresponds to the element at free end of the workpiece when it extends to infinity, i.e., the element
80
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
Fig. 6. The wave generated at different locations of the workpiece after throw-off element consideration.
Location number
1
2
3
4
Distance away from the irradiated spot center (m)
0.125
0.25
0.375
0.5
behaves as a throw-off solid piece in that it is conduit for energy out of the system. Consequently, no reflection of the wave from the free end of the workpiece is encountered. Fig. 6 shows the wave amplitude with time after throw-off element consideration. In the first location, the wave has a certain pattern after 4 ms: In the second location, this pattern appears after 6 ms; and in the next location the regular pattern is replaced by the decaying wave pattern. Consequently, despite the dispersion effect of the workpiece material, the flexural wave generated has a certain pattern at locations close to the irradiated spot center. This identity disappears as the location moves away from the irradiated spot center. The peak amplitude of the
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
81
wave at different locations varies, i.e. the peak amplitude at location 2 is higher than that corresponding to the first location. This may be because of the partial overlapping of the second mode with the first mode, which reduces the peak amplitude of the wave. It should be noted that two modes of wave propagation occur. The second mode is faster than the first mode, so it is always ahead of the first mode. As the wave propagates in the workpiece towards the free end, the two modes are expected to separate because of the difference in their speed. However, the distance along the workpiece is not large enough to distinguish both modes of the travelling waves at present.
4. Conclusions The thermal and elastic response of steel due to laser heating pulse is investigated. The normal component of the thermal stress is considered as the source for the flexural wave generation in the workpiece material. The travelling flexural wave is computed at four locations at the workpiece surface. In general, the normal component of the thermal stress in the surface vicinity is tensile and no unique pattern of the flexural wave is obtained at different locations on the surface of the workpiece. This is due to the dispersion effect of the workpiece material, overlapping of two modes, and interference of the reflected waves emanating from the free ends of the workpiece. The specific conclusions derived from the present study can be listed as follows: 1. The temperature rises rapidly at the surface in the early heating period and the rate of temperature rise reduces as the heating progresses. The temperature decay rate in the beginning of the cooling cycle is higher than that corresponding to high cooling periods. This is because the temperature gradient in the surface vicinity of the workpiece reduces at a slow rate as the cooling period progresses. 2. The normal component of the thermal stress is tensile in the surface vicinity and it varies with time such that two peaks are observed, one in the heating cycle and the other in the cooling cycle. The strain continues to increase in the early cooling period provided that it reduces rapidly as the cooling period progresses further. 3. The amplitude of the travelling wave reduces as the location at the workpiece surface moves towards the free end of the workpiece. This is because of the damping and dispersion effects of the workpiece material and interference of the reflected waves from the free ends of the workpiece. The pattern of the travelling wave changes towards the workpiece end and the wave amplitude breaks into high frequency peaks. This is because of the reflected waves and the partial overlapping of the two modes of the travelling wave. 4. The throw-off analysis indicates that the reflected waves have significant effect on the travelling wave, since the time required for the reflected wave to reach the center of the irradiated spot is in the order of 0:16 10 3 s; which is considerably less than the time domain considered in the present study. Moreover, a certain pattern is identified in the travelling wave at locations close to the center of the
82
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
irradiated spot when throw-off element is introduced at the free end of the workpiece. This pattern is shifted in time domain as the location moves away from the irradiated spot center. This is because of the effect of the second mode on the travelling wave; in which case, the second mode travels faster than the first mode and the partial overlapping of two modes occurs, which modifies the wave pattern.
Acknowledgements The authors acknowledge the support of King Fahd University of Petroleum and Minerals, Dhahran, Saudi Arabia, for this work.
References [1] Yilbas BS, Sami M, Shuja SZ. Laser-induced thermal stresses on steel surface. Opt Lasers Eng 1998;30:25–37. [2] Wang HG, Guan YH, Chen TL, Zhang JT. A study of thermal stresses during laser quenching. J Mater Proc Technol 1997;63:550–3. [3] Li K, Sheng P. Plane stress model for fracture of ceramics during laser cutting. Int J Mach Tools Manuf 1995;35:1493–506. [4] Elperin T, Rudin G. Thermo-elasticity problem for a multilayer coating-substrate assembly irradiated by a laser beam. Int Commun Heat Mass Transfer 1996;23:133–42. [5] Modest MF. Transient elastic and viscoelastic thermal stresses during laser drilling of ceramics. ASME J Heat Transfer 1998;120:892–8. [6] Schulz D, Becker D, Franke J, Kemmerling R, Herziger G. Heat conduction losses in laser cutting of metals. J Phys D 1993;26:1357–63. [7] Shuja SZ, Yilbas BS. Pulsative heating of surfaces. Int J Heat Mass Transfer 1998;41:3899–918. [8] Neto D, Lima C. Non-linear three-dimensional temperature profiles in pulsed laser heated solids. J Phys D 1994;27:1795–804. [9] Yilbas BS. Laser heating process and experimental validation. Int J Heat Mass Transfer 1997;40(5):1131–43. [10] Mason WP, Thurston RN. Physical acoustics, vol. 18. San Diego: Academic Press, 1988. p. 21–123. [11] Scruby CB, Drain LE. Laser ultrasonics; techniques and applications. Bristol, UK: Adam Hilger, 1990. [12] Hsu NN. US National Bureau of Standarts Report, NBSIR85-3234, 1985. [13] Auld BA. Acoustic fields and waves in solids, vol. 1. Florida: Krieger, 1990. [14] Dubois M, Enguehard F, Bertland L. Modeling of laser thermoelastic generation of ultrasound in an orthotropic medium. Appl Phys Lett 1994;65:554–6. [15] Cheng J, Wu L, Zong S. Thermoelastic response of pulsed photothermal deformation of thin plates. J Appl Phys 1994;76:716–22. [16] Graff KF. Wave motion in elastic solids. New York: Oxford University Press, 1975. [17] Dewhurst RJ, Edwards C, Mckie ADW, Palmer SB. Estimation of thin metal sheet using laser generated ultrasound. Appl Phys Lett 1987;51(14):1066–8. [18] Cheng J, Bertelot Y. Theory of laser generated transient Lamb waves in orthotropic paltes. J Phys D 1996;29:1857–67. [19] Shuja SZ, Yilbas BS. Gas-assisted repetitive pulsed heating of a steel surface. Proc Inst Mech Eng, Part C 1998;212:741–57. [20] Gotski TB, Hussaini MY, Lumley JL. Simulation and modeling of turbulent flows. New York: Oxford University Press, 1996.
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
83
[21] Incropera FP, Dewitt DP. Introduction to heat transfer F appendix. New York: Wiley, 1985. p. 667–96. [22] Versteeg HK, Malalasekera W. An introduction to computational fluid dynamics: the finite volume method. New York: Longman Scientific and Technical, 1995. [23] Shuja SZ. Simulation of three-dimensional gas-assisted heating of solid substance. Ph.D. thesis, Mechanical Engineering Department, KFUPM, Dhahran, Saudi Arabia. [24] Patankar SV. Computer analysis of fluid flow and heat transfer. Swansea: Pinridge, 1981. p. 223–52 [Chapter 8]. [25] Bathe KJ. Finite element procedures. Englewood Cliffs, NJ: Prentice-Hall, 1996. p. 148–92 [Chapter 6]. [26] Doyle JF. Wave propagation in structures. New York: Springer, 1989. [27] Mahmood F. Transient response of varying thickness Timonshenko beam under poind load. M.Sc. thesis, Mechanical Engineering Department, KFUPM, Dhahran, Saudi Arabia.
Laser pulse heating of steel surface and flexural wave analysis B.S. Yilbas*, M. Faisal, S.Z. Shuja, A.F.M. Arif King Fahd University of Petroleum and Minerals, Dhahran, Saudi Arabia Received 13 July 2000; received in revised form 1 June 2001; accepted 1 August 2001
Abstract Laser gas-assisted material processing finds wide application in industry. The modelling of heating, elastic response of the substrate material, and the wave analysis gives insight into the laser workpiece interaction. In the present study, laser gas-assisted heating of steel is considered. The normal component of the thermal stress is taken as the source of load for the flexural wave generation in the material. The flexural wave generated is simulated and the wave characteristics are analyzed at four locations at the workpiece surface. The numerical scheme employing a control volume approach is introduced when solving the governing equations of flow and heat transfer while finite element and spectran element methods are used when solving the stress and wave equations. It is found that the normal component of the stress is tensile. The dispersion effect of the workpiece material, interference of the reflected beam, and partial overlapping of second mode of the travelling wave enable to identify a unique pattern in the travelling wave in the substrate. r 2001 Elsevier Science Ltd. All rights reserved. Keyword: Laser heating flexural wave generation
1. Introduction Lasers are widely used in industry due to their precision of operation, high speed of processing, and low cost. When the laser beam is absorbed by the substrate material, temperature of the irradiated region rises, which in turn results in thermal stress generation due to thermal expansion of the workpiece material. Thermal loading produces flexural waves propagating in the substrate. The workpiece acts as *Corresponding author. Tel.: +966-3-860-2540; fax: +966-3-860-2949. E-mail address: [email protected] (B.S. Yilbas). 0143-8166/02/$ - see front matter r 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 1 4 3 - 8 1 6 6 ( 0 1 ) 0 0 1 4 9 - X
64
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
Nomenclature a Gaussian parameter A cross-sectional area of workpiece, m2 cp specific heat capacity, J=kg K C1 ; C2 ; Cm coefficients in the k–e turbulence model d diameter of the nozzle, m D elasticity matrix Dp surface displacement, m f1 ; f2 ; fm coefficients in the low Reynolds number, k–e model Df frequency step, Hz E modulus of elasticity, MPa G rate of generation of k; W=m3 I0 peak power intensity, W=m2 I moment of inertia of cross section, m4 k turbulent kinetic energy, W=m3 k wave number (Eqs. 26–28) M bending moment, Nm N number of data points in power spectrum p pressure, Pa Pr varible Prandtl number (function of temperature) q heat flux, W=m2 Re Reynolds number r distance in the radial direction, m S unsteady spatially varying source (Eq. (9)), W=m3 t time, s T temperature, K U arbitrary velocity, m=s V radial velocity, m=s W axial velocity, m=s x arbitrary direction, m zn distance to the nearest wall, m z distance in the axial direction, m DTðr; z; tÞ Temperature rise at a point ðr; zÞ at time¼ t with respect to that at t ¼ 0 Tref Reference temperature at t ¼ 0 ff B g the applied body force fPg the applied pressure vector % fFg concentrated nodal forces to the element fdUg virtual displacement fdUs g virtual displacement on the boundary where pressure is prescribed % fdUg virtual displacement of boundary nodes where concentrated load is prescribed ½B strain displacement gradient matrix % fUg nodal displacement vector
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
65
½N matrix of shape (or interpolation) functions R fF th g element thermal load vector, 8 ½BTR½D½eth d8 T B fFgb element applied body force 8 ½N ½f d8 R vector, s T fFg element pressure vector, R ½Nn ½P d ½Ke element stiffness matrix, 8 ½BT ½D½B d8 Greek letters s feg feth g e0 a u G d e Z l k k m
r F f f o
stress vector, MPa total strain vector thermal strain vector equivalent strain coefficient of thermal expansion, 1/K Poisson’s ratio diffusion coefficient, m2 =s absorption coefficient, 1=m energy dissipation, W=kg damping coefficient thermal conductivity, W=mK turbulence intensity shear correction factor (Eqs. 22; 23) variable dynamic viscosity (function of temperature for laminar), N s=m2 effective viscosity, ðml þ mt Þ; N s=m2 variable kinematic viscosity (function of temperature for laminar), m2 =s density (function of temperature and pressure for gas), kg=m3 viscous dissipation, W=m3 arbitrary variable bending slope circular frequency, Hz
Subscript amb i; j jet l max t
ambient arbitrary direction conditions at jet inlet laminar maximum turbulent
me n
a dispersive media influencing the wave propagation. The wave response of the workpiece is governed by three system parameters namely, damping, mass, and stiffness. Consequently, the geometry, the mechanical properties, and the position of the workpiece determine the wave propagation characteristics.
66
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
In order to analyze the laser-induced flexural waves, the stress field generated in the laser-irradiated region needs to be analyzed first. This is because the load generating the flexural wave is due to the axial and transverse stress components in the axisymmetric system without torsional motion. Considerable research studies were conducted on the laser-induced thermal stresses in solids. The laser-induced thermal stresses on steel surfaces were investigated by Yilbas et al. [1]. They showed that excessive stress levels occurred in the irradiated spot, which in turn resulted in the deformation at the surface of the substrate material. The thermal stresses generated during the laser quenching was investigated by Wang et al. [2]. They indicated that the residual stress developed in the surface vicinity resulted in the hardening of the surface. The experimental and numerical studies of fracture initiation in thin aluminum oxide ceramics during laser cutting were carried out by Li and Sheng [3]. They introduced plane stress model when predicting the stresses along the cut edges. They showed that low cutting speed reduces the stress level and minimizes the possible fracture initiation at the cut edges. An analytical solution of a two-dimensional thin coating attached to a substrate was obtained by Elperin and Rudin [4]. They indicated that the temperature difference across the coating generated excessive stress levels, which in turn caused inelastic deformation of the coating. Modest [5] investigated the elastic and viscoelastic stresses during laser drilling of ceramics. He indicated that the material softened in the region closer to the melt zone. The temperature field generated in the laser-irradiated region plays an important role in stress development and flexural wave initiation in the substrate. Considerable research studies were carried out to investigate the temperature field in the substrate during laser heating process. Heat conduction losses in laser cutting of metals was investigated by Schultz et al. [6]. They showed that the theoretical predictions agreed well with the experimental findings. Shuja and Yilbas [7] investigated the gas-assisted laser heating of steel substrate. They introduced the three-dimensional conjugate heating model and predicted the temperature field in both the solid and the gas side. They indicated that the impinging gas jet had little influence on the surface temperature rise. Neto Diniz and Lima [8] employed a three-dimensional parabolic heat diffusion equation to predict the temperature field during laser pulse heating of steel surfaces. They indicated that the oscillation in the pulse had no significant effect on the resulting heat-affected zone depth. Yilbas [9] introduced a three-dimensional heating model with a stationary Gaussian laser beam source. The phase change process was employed in the analysis. He indicated that the temperature below the surface became higher than that of the surface as the rate of surface evaporation increased. The laser-induced waves in solids received considerable attention, since it served as a nondestructive tool for materials evaluation [10,11]. Several analytical methods were introduced when modelling the thermoelastic process that took place during the laser heating pulse. Some of these include Green function method, and Laplace and Henkel transformations [12,13]. However, the thermal stresses developed in the substrate vary with time and it is, in general, difficult to present in a simple mathematical function. Therefore, the analytical solution to the transient problem
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
67
with nonordered input load becomes very difficult to obtain. Consequently, the general trend is to solve the governing equations numerically. Dubois et al. [14] investigated numerically the laser thermoelastic generation of ultrasound in an orthotropic medium. They compared the predictions with the experimental results and indicated that their predictions agreed well with the experimental findings. They further argued that the model developed could be used to optimize the laser generation of ultrasound in various materials. Cheng et al. [15] studied the dynamic thermoelastic response of pulsed photothermal deformation deflection detection for Q-switched laser pulses. They indicated that when the laser pulse rise time decreases in the order of 10 ns or less, the pulsed photothermal deformation deflection signal reflected a totally dynamic wave behavior. The laser-induced waves such as Lamb waves can also be used to determine the properties of the substrate material. The Lamb wave is generated due to normal and shear loads applied to the surface of the half-space [16]. Dewhurst et al. [17] introduced a Lamb wave technique to describe the measurement of the film thickness on the substrate material. A theoretical study for modelling the thermoelastic excitation of transient Lamb waves propagating along the principal direction in an orthotropic plate was carried out by Cheng and Berthelod [18]. They provided a quantitative analysis for noncontact and nondestructive detection of the elastic stiffness properties of the machine-made paper by using the laser-generated Lamb wave technique. In the light of the above arguments, the present study is carried out to investigate the flexural wave propagation in the substrate during the laser heating process. The thermal stresses developed in the substrate are predicted first and the wave generated due to axial stress component in the surface vicinity of the substrate is analyzed using a spectral finite element method. Since the workpiece has a certain length, the reflecting waves from the free supported ends of the workpiece are also considered in the analysis.
2. The mathematical model The gas-assisted laser pulse heating process is considered when modelling the laser workpiece interaction provided that conduction limited heating case is taken into account in the present study [19]. Fig. 1 shows the schematic view of the heating process. 2.1. Flow and heat conduction equations The unsteady flow equations need to be solved to obtain the flow field due to axisymmetric gas jet impingement. The continuity and momentum equations are qr q þ ðrUi Þ ¼ 0 qt qxi
ð1Þ
68
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
Nozzle
r Nozzle exit 1 mm dia. Air gas jet velocity100 m/s Stand off 2 mm
z
Solid gas interface
Solid thickness 2 mm Stress free boundary Free support
Stress free boundary Free support
Fig. 1. A view of impinging gas and laser workpiece arrangement.
and
q q qp q qUj ðrUj Þ þ ðrUi Uj Þ ¼ þ ðm þ mÞ ; qt qxi qxj qxi t qxi
ð2Þ
where mt is the eddy viscosity which has to be specified by a turbulence model. In the present study, low Reynolds number k–e turbulence model is introduced to account for the turbulence effect of the impinging gas, which is air [19]. The partial differential equation governing the transport of thermal energy has the form qT q q mt m qT þ ðUi TÞ ¼ þ : ð3Þ qt qxi qxi Prt Pr qxi
2.1.1. Boundary conditions for flow equations Laminar boundary conditions are set for the mean-flow variables, and the boundary conditions k ¼ 0 and de=dz ¼ 0 are applied at the wall. The low-Reynolds number extension does not employ wall functions; therefore, the grid employed normal to the main flow direction needs to be distributed so as to give a high concentration of grid cells near the wall, with the wall-adjacent node positioned at zþ ¼ rzun =mp1:0: Inlet to control volume. Ui ¼ specified
and
T ¼ constant:
ð4Þ
The kinetic energy of turbulence is estimated according to some fraction of the square of the average inlet velocity [20] k ¼ lu%2 ; where u% is the average inlet velocity and l is a fraction.
ð5Þ
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
69
The dissipation is calculated according to the equation [11] k3=2 ; ð6Þ bd where d is the inlet diameter. The values l ¼ 0:03 and b ¼ 0:005 are commonly used and may vary slightly in the literature [20]. Outlet to control volume It is assumed that the flow extends over a sufficiently long domain; therefore, it is fully developed at the exit section. Thus, for any variable f the condition is qðrfÞ ¼ 0; ð7Þ qx where x is the arbitrary outlet direction. Symmetry axis The radial derivative of the variables is set to zero at the symmetry axis, i.e. qf ¼ 0 and V ¼ 0: ð8Þ qr Solid fluid interface The temperature at the solid–gas interface is considered as the same, i.e. qTwgas qTwsolid Twsolid ¼ Twgas and Ksolid ¼ Kgas : qz qz The unsteady heat conduction equation. q q qT ðcp rTÞ ¼ K þ S; ð9Þ qt qxi qxi e ¼ Cm
where S is the unsteady spatially varying laser output power intensity distribution and is considered as Gaussian with 1=e2 points equal to 0:375 mm; from the center of the beam. Therefore, S is stated as follows: 2 I0 r S ¼ pffiffiffiffiffiffiffiffi exp 2 d expð dzÞf ðtÞ; ð10Þ a 2pa pffiffiffiffiffiffiffiffi where ðI0 = 2paÞ expðr2 =a2 Þ is the intensity distribution across the surface, expð dzÞ is the absorption function, and f ðtÞ is the function accommodating the time variation of the pulse shape. The pulse properties are given in Table 1. Boundary conditions for the heat conduction equation Convection with a constant coefficient for still air is considered at the z ¼ zmax boundary for the plate. The continuity of temperature between the solid and the gas is enforced at the interface;
Table 1 Pulse properties used in the simulation Rise time (ms)
Fall time (ms)
Total pulse length (ms)
Peak power intensity ðW=m2 Þ
0.2
1.08
1.48
109
70
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
and a constant temperature, T ¼ Tamb ; is assumed for the distance far away from the laser source. 2.1.2. Variable properties An equation of state is used for the impinging gas and the specific heat capacity and thermal conductivity for both air and steel were considered only as a function of temperature. The temperature dependence of properties is given in [21]. 2.1.3. Calculation of the flow field variables The control volume approach is used when solving the governing equations of flow and heat transfer numerically. The differential equation is integrated over the control volume to yield the discretization equation. The main reasons for choosing the control-volume formulation for flow field are its simplicity and easy physical interpretation [22]. The discretization process is not given here due to lengthy arguments, but refer to [23]. The grid used in the present calculations has 38 70 mesh points, provided that the grid independent test is satisfied. The details of the grid orientation and grid independent test results are given in [9]. The two problems of determining the pressure and satisfying continuity are overcome by adjusting the pressure field so as to satisfy continuity. This arrangement gives a convenient way of handling the pressure linkages through the continuity equation and is known as the SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algorithm. The details of this algorithm are found in [24]. The governing equation for heat conduction in solid (Eq. (9)) can be written in the form of flow equations. Thus, the discretization procedure leads to algebraic equations of the form similar to flow equations with temperature T replacing the general variable. 2.2. Thermal stress modelling The temperature field results in the thermal stresses being generated in the substrate, which can lead to the elastic–plastic displacement in the substrate material. The stress is related to strains by fsg ¼ ½Dfee g;
ð11Þ
where fsg is the stress vector, and ½D is the elasticity matrix. fee g ¼ feg feth g; where feg is the total strain vector and feth g is the thermal strain vector. Eq. (11) may also be written as feg ¼ ½D 1 fsg þ feth g:
ð12Þ
Since the present case is axially symmetric, and the material is assumed to be isotropic, the above stress–strain relations can be written in cylindrical
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
coordinates as 1 err ¼ ½srr uðsyy þ szz Þ þ aDTðr; z; tÞ; E 1 eyy ¼ ½syy uðsrr þ szz Þ þ aDTðr; z; tÞ; E 1 ezz ¼ ½szz uðsrr þ syy Þ þ aDTðr; z; tÞ; E 1 erz ¼ srz : G
71
ð13Þ
where E; u; and a are the modulus of elasticity, Poisson’s ratio, and coefficient of thermal expansion, respectively. DTðr; z; tÞ represents the temperature rise at a point ðr; zÞ at time ¼ t with respect to that at t ¼ 0 corresponding to a stress-free condition. A typical component of thermal strain from Eq. (13) is eth ¼ aDTðr; z; tÞ ¼ aðTðr; z; tÞ Tref Þ; where Tref is the reference temperature at t ¼ 0: When a is a function of temperature, then Eq. (14) becomes Z T th e ¼ aðTÞ dT:
ð14Þ
ð15Þ
Tref
The present study uses a mean or weighted-average value of a; such that eth ¼ a% ðTÞðT Tref Þ;
ð16Þ
where a% ðTÞ is the mean value of coefficient of thermal expansion and is given by RT Tref aðTÞ dT a% ðTÞ ¼ : ð17Þ Tðr; z; tÞ Tref To develop a finite element procedure for stress computation, the standard displacement-based finite element method is used. The basis of this approach is the principle of virtual work, which states that the equilibrium of any body under loading requires that for any compatible small virtual displacements (which are zero at the boundary points and surfaces and correspond to the components of displacements that are prescribed at those points and surfaces) imposed on the body in its state of equilibrium, the total internal virtual work or strain energy ðdUÞ is equal to the total external work due to the applied thermally induced loads ðdVÞ; i.e. dU ¼ dV: For the static analysis of problems having linear geometry and thermoelastic material behavior, one can derive the following equation using standard procedure [25]. Z ðfdegT ½Dfeg fdegT ½Dfeth gÞ d8 8 Z Z X % T fFg; % ¼ fdUgT ff B g d8 þ fdUs gT fPg dO þ fdUg ð18Þ 8
72
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
% the where ff B g is the applied body force, fPg the applied pressure vector, fFg concentrated nodal forces to the element, fdUg the virtual displacement, fdUs g the % the virtual displacement on the boundary where pressure is prescribed, and fdUg virtual displacement of boundary nodes where concentrated load is prescribed. The strains may be related to the nodal displacement by % feg ¼ ½BfUg;
ð19Þ
% the nodal displacement where ½B is the strain displacement gradient matrix, and fUg vector. The displacements within the elements are related to the nodal displacement by % fUg ¼ ½NfUg; ð20Þ where ½N is the matrix of shape (or interpolation) functions. Eq. (18) can be reduced to the following matrix form: % fF th g ¼ fFgb þ fFgs þ fFg % ½Ke fUg ð21Þ R R T T where ½Ke ¼ 8 ½B ½D½B d8 is the element stiffness fF th g ¼ 8 ½B ½D R matrix, b T B th ½A d8 the element Rthermal load vector, fFg ¼ 8 ½N ½f d8 the element body force vector, fFgs ¼ ½Nn T ½P d the element pressure vector, and ½Nn ¼ matrix of shape functions for normal displacement at the boundary surface. Assembly of element matrices and vectors of Eq. (21) yields % ¼ fRg; % ½Kfdg % and fRg % are the global stiffness matrix, global nodal displacement where ½K; fdg vector, and global nodal load vector, respectively. Solution of the above set of simultaneous algebraic equations gives unknown nodal displacements and reaction forces. 2.3. Wave analysis In the analysis, the workpiece is divided into a number of elements and each element can be treated as uniform. The displacement vector is then transformed from time domain into frequency domain using fast Fourier transformation (FFT). The transform response is fed back into the governing differential equations as solved at each frequency component using finite element method (FEM) and the transient response is recovered by using the inverse FFT. 2.3.1. Formulation The solution in the wave problems is general function of space and time. If the time variation of the solution is focussed on a particular point in space, then it can be presented by a spectral function. The spectral formulation starts with the equations of motion of the workpiece (as free supported beam) including inertia terms. After dividing the workpiece into a number of elements, cross sectional area and second moment of area in the center of each element can be considered as constant
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
73
over the entire element. Therefore, the equation of motion can be written as [26] 2 q y qf q2 y GAk ð22Þ ¼ rA 2 q 2 qx qx qt and q2 f qy q2 f f ¼ rI þ GAk : qx2 qx qt2
ð23Þ
The displacement and the bending slope can be written in the spectral form as [26] X ; yðx; tÞ ¼ yn ðx; on Þeion t ð24Þ n
and fðx; tÞ ¼
X
;
fn ðx; on Þeion t :
ð25Þ
n
The spectral components yn ; and fn ; have the solutions ;
yn ¼ A1 e ik1 x þ B1 e ik2 x þ C1 eik1 x þ D1 eik2 x
ð26Þ
and ;
fn ¼ A2 e ik1 x þ B2 e ik2 x þ C2 eik1 x þ D2 eik2 x ;
ð27Þ
where k1 and k2 are the wave numbers and A1 ; A2 ; B1 ; B2 ; C1 ; C2 and D1 ; D2 are the frequency-dependent coefficients, which are complex in nature. The first term in Eqs. (26) and (27) represents the waves moving in the forward direction while the last two describe the backward moving waves. Substituting yn ; and fn ; into Eqs. (22) and (23), and considering the load to be applied at the modes, yields GAkð k21 A1 þ ik1 A2 Þ ¼ o2 rAð A1 Þ; GAkð k21 B1 þ ik1 B2 Þ ¼ o2 rAð B1 Þ; GAkð k21 C1 þ ik1 C2 Þ ¼ o2 rAð C1 Þ; GAkð k21 D1 þ ik1 D2 Þ ¼ o2 rAð D1 Þ; where the coefficients are as follows:
A2 ¼
o2 rA þ k21 GAk A1 ; ik1 GAk
B2 ¼
o2 rA þ k21 GAk B1 ; ik1 GAk
ð28Þ
74
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
2 o rA k21 GAk C1 ; ik1 GAk 2 o rA k21 GAk D2 ¼ D1 : ik1 GAk
C2 ¼
The shear force and momentum equations can be written as 2
;
d f Q’ ¼ EI 2 o2 rI f dx
ð29Þ
and ;
’ ¼ EI d f : M dx
ð30Þ
The dynamic stiffness is obtained by first relating the coefficients to the nodal displacement as ;
; ;
;
;
½A; B; C; D ¼ a ½y1 ; f1 ; y2 ; f2 ;
ð31Þ
where ½a; ¼ ½b and ½b 1 are given in [27]. The nodal loads are obtained by using the shear force and momentum equations Eqs. (29) and (30), i.e. ;
;
½F T ¼ ½k ½dT ; where ;
;
;
;
;
½F T ¼ ½Q1 M1 y Qn Mn and ;
;
;
;
½dT ¼ ½y1 f1 y y1 f1 where n is the number of nodes. 2.3.2. Initial and boundary conditions Initially the displacements and velocities are set to zero. The free supported boundary condition at the workpiece ends is considered. In this case, at free ends of the workpiece ðx ¼ 0 ¼ 1Þ; the bending moment is set to zero at all frequencies, i.e. d f; ðI; oÞ ¼ 0: dx In order to observe the effect of reflecting waves from the free ends of the workpiece, one end of the workpiece is assumed to be at infinity. This boundary condition (element extends to infinity) requires that the element behaves as a throwoff beam, which acts as a conduit for energy out of the system.
75
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
2.3.3. Method of solution In order to obtain the transient response of the workpiece using a spectral finite element method, the input excitation signal, which is the load generated due to normal component of the thermal stress, at discrete time steps needs to be converted into frequency domain by FFT. The sampling rate Dt must be in accordance with Nyquist frequency of the signal. The frequency step is obtained from 1 Df ¼ ; NDt where N is the number of data points in the power spectrum. The signal is bondlimited between the minimum frequency ðfmin Þ and the maximum frequency ðfmax Þ: Beginning from the minimum frequency, the assembled dynamic stiffness matrix and the assembled load vector in frequency domain are obtained by the finite element method, i.e.: ;
; ;
;
;
½F ¼ ½K ½d
ð32Þ
or ;
½d ¼ ½K 1 ½F : After applying the boundary conditions to Eq. (32), ½K is inverted by Gauss reduction method. Consequently, the response is obtained at fmin : This procedure is repeated for other frequencies upto fmax : Hence, the displacement vector ½d ; is obtained at the desired element. The frequency-dependent response is reconstructed with time domain by using inverse Fourier transformation method. In order to improve the stability of the computation, damping is introduced through the wave number k in the form k ¼ k0 ð1 iZÞ where k0 is the undamped value and Z is taken as 0.025. The material properties are given in Table 2 while workpiece size and number of elements used in the simulation are given in Table 3. Table 2 Properties of steel used in the simulation and the number of elements employed E (GPa)
G (GPa)
r ðkg=m3 Þ
k
u
Z
207
77.6
7836
0.67
0.3
0.025
Table 3 Workpiece size and number of elements used in simulations Length (m)
Width (m)
Thickness (m)
Number of elements
1
0.2
0.002
8
76
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
3. Results and discussions The thermal stresses inside the substrate due to the temperature field developed by a laser heating pulse is modelled. The flexural wave generation in the substrate material due to the normal component of the thermal stress is analyzed. In order to investigate the wave propagation along the transverse direction in the workpiece, six locations with three located at each side from the center of the workpiece along the transverse direction are considered. This is necessary, since the workpiece material disperses the travelling waves in amplitude and frequency domain. The wave reflections from both ends of the workpiece also modify the travelling wave properties. Consequently, to analyze the influence of reflected waves emanating from the ends of the workpiece, a case in which one end of the workpiece is located at infinity is considered. This provides no reflection of the wave from one end of the workpiece. Fig. 2 shows the temporal variation of the surface temperature at the center of the heated spot (at r ¼ 0) for convective and nonconvective boundary conditions at the surface. The temperature rises rapidly in the early heating period and the rate of rise of surface temperature reduces as the heating progresses. In the cooling cycle, the decay rate of the surface temperature is faster in the beginning of the cooling cycle and the decay rate slows as the cooling period progresses. Consequently, the material response to the heating pulse changes in the cooling cycle as compared to that in the heating cycle. This is because once the laser source disappears, the rate of energy diffused from the surface vicinity of the substrate to the bulk through the conduction becomes large immediately after the laser pulse ends, i.e. the large temperature gradient attainment in the surface region during the heating cycle accelerates the conduction cooling of the surface vicinity. In the case of the long cooling period, the
Fig. 2. Temporal variation of surface temperature at the center of the irradiated spot with and without convective boundary condition at the surface.
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
77
temperature gradient in the surface vicinity is low and the energy diffusion from the surface vicinity to the bulk of the substrate slows down. When comparing temperature profiles with and without convective boundary condition at the surface, it is evident that both temperature profiles are identical. This indicates that convective cooling of the surface is negligibly smaller than the energy input to the substrate material due to the absorption of a laser beam as consistent with the previous work [9]. Fig. 3a shows the temporal variation of the normal (axial) component of the thermal stress while Fig. 3b shows the temporal variation of the equivalent strain in the surface vicinity. The equivalent strain increases as the heating progresses
25
STRESS (MPa)
18
11
4
-3
-10 0
0.5
1
1.5
2
1.2
1.6
TIME (ms)
(a) 0.025
STRAIN
0.020
0.015
0.010
0.005
0.000 0
(b)
0.4
0.8
TIME (ms)
Fig. 3. (a) Temporal variation of load generated in the surface vicinity of the workpiece due to normal component of the thermal stress; (b) Temporal variation of strain in the surface vicinity of the workpiece.
78
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
provided that the increase in strain is continuous even in the cooling cycle. It reduces rapidly as the cooling period progresses further. This is because the elastic response of the material is not as fast as the temperature field, i.e. the expansion is continuous for a while after the laser source vanishes. The normal component of the stress varies considerably with time. It results in two maxima; one in the heating cycle and the other in the cooling cycle. The normal component of the stress in the surface vicinity is tensile and it reaches as high as 20 MPa; which is less than the yield stress of the workpiece material. Fig. 4a shows the load variation with time while Fig. 4b shows the power spectrum density curve. The temporal variation of the load is similar to the normal stress component as shown in Fig. 3a, i.e. the value of the stress component is multiplied by the size of the irradiated spot to obtain the load. Since the total power in a signal is constant, the low and high frequency limits of the signal are determined from Fig. 4b. In the present case, Nyquist frequency (fc ¼ 1=2Dt; where Dt is the sampling rate) is selected as 2000. It should be noted that fc is the maximum frequency contained in the signal. Fig. 5 shows the flexural wave amplitude variation with time at four equally spaced locations at the workpiece surface. The first location is 12:5 cm away from the center of the heated spot and following locations are 12:5 cm away from the first location, respectively. The wave amplitude varies with time and no dominant pattern is observed in the wave. This is because the initial wave form at the center of the workpiece is modified as it travels along the workpiece due to: (i) the dispersion effect of the workpiece material, (ii) overlapping of wave modes, and (iii) the reflected waves from the free ends of the workpiece that interfere with the travelling wave. The group speed of the first wave mode is about 3125 m=s; therefore, the reflection from the free ends takes about 0:16 10 3 s to reach the center of the workpiece. Consequently, the reflected wave modifies the travelling wave in the time and amplitude domain. As the location from the first point changes to the following
Fig. 4. (a) Temporal variation of load generated due to normal component of the thermal stress; (b) Power spectrum density of thermal load.
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
79
Fig. 5. The wave generated at different locations of the workpiece.
Location number
1
2
3
4
Distance away from the irradiated spot center (m)
0.125
0.25
0.375
0.5
point along the surface, the amplitude of the wave decreases and the travelling wave pattern changes. This indicates that the travelling wave vanishes as it travels towards the free ends of the workpiece, i.e. in the present study, damping coefficient of the workpiece material is taken as 0.025. Moreover, at locations close to the free ends of the workpiece, the wave appears as in high frequency mode. This is because the travelling wave has a low amplitude in these locations and reflected waves modify easily the travelling wave pattern. In order to identify the influence of the reflected waves on the travelling wave, throw-off element is introduced at one end of the workpiece. This corresponds to the element at free end of the workpiece when it extends to infinity, i.e., the element
80
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
Fig. 6. The wave generated at different locations of the workpiece after throw-off element consideration.
Location number
1
2
3
4
Distance away from the irradiated spot center (m)
0.125
0.25
0.375
0.5
behaves as a throw-off solid piece in that it is conduit for energy out of the system. Consequently, no reflection of the wave from the free end of the workpiece is encountered. Fig. 6 shows the wave amplitude with time after throw-off element consideration. In the first location, the wave has a certain pattern after 4 ms: In the second location, this pattern appears after 6 ms; and in the next location the regular pattern is replaced by the decaying wave pattern. Consequently, despite the dispersion effect of the workpiece material, the flexural wave generated has a certain pattern at locations close to the irradiated spot center. This identity disappears as the location moves away from the irradiated spot center. The peak amplitude of the
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
81
wave at different locations varies, i.e. the peak amplitude at location 2 is higher than that corresponding to the first location. This may be because of the partial overlapping of the second mode with the first mode, which reduces the peak amplitude of the wave. It should be noted that two modes of wave propagation occur. The second mode is faster than the first mode, so it is always ahead of the first mode. As the wave propagates in the workpiece towards the free end, the two modes are expected to separate because of the difference in their speed. However, the distance along the workpiece is not large enough to distinguish both modes of the travelling waves at present.
4. Conclusions The thermal and elastic response of steel due to laser heating pulse is investigated. The normal component of the thermal stress is considered as the source for the flexural wave generation in the workpiece material. The travelling flexural wave is computed at four locations at the workpiece surface. In general, the normal component of the thermal stress in the surface vicinity is tensile and no unique pattern of the flexural wave is obtained at different locations on the surface of the workpiece. This is due to the dispersion effect of the workpiece material, overlapping of two modes, and interference of the reflected waves emanating from the free ends of the workpiece. The specific conclusions derived from the present study can be listed as follows: 1. The temperature rises rapidly at the surface in the early heating period and the rate of temperature rise reduces as the heating progresses. The temperature decay rate in the beginning of the cooling cycle is higher than that corresponding to high cooling periods. This is because the temperature gradient in the surface vicinity of the workpiece reduces at a slow rate as the cooling period progresses. 2. The normal component of the thermal stress is tensile in the surface vicinity and it varies with time such that two peaks are observed, one in the heating cycle and the other in the cooling cycle. The strain continues to increase in the early cooling period provided that it reduces rapidly as the cooling period progresses further. 3. The amplitude of the travelling wave reduces as the location at the workpiece surface moves towards the free end of the workpiece. This is because of the damping and dispersion effects of the workpiece material and interference of the reflected waves from the free ends of the workpiece. The pattern of the travelling wave changes towards the workpiece end and the wave amplitude breaks into high frequency peaks. This is because of the reflected waves and the partial overlapping of the two modes of the travelling wave. 4. The throw-off analysis indicates that the reflected waves have significant effect on the travelling wave, since the time required for the reflected wave to reach the center of the irradiated spot is in the order of 0:16 10 3 s; which is considerably less than the time domain considered in the present study. Moreover, a certain pattern is identified in the travelling wave at locations close to the center of the
82
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
irradiated spot when throw-off element is introduced at the free end of the workpiece. This pattern is shifted in time domain as the location moves away from the irradiated spot center. This is because of the effect of the second mode on the travelling wave; in which case, the second mode travels faster than the first mode and the partial overlapping of two modes occurs, which modifies the wave pattern.
Acknowledgements The authors acknowledge the support of King Fahd University of Petroleum and Minerals, Dhahran, Saudi Arabia, for this work.
References [1] Yilbas BS, Sami M, Shuja SZ. Laser-induced thermal stresses on steel surface. Opt Lasers Eng 1998;30:25–37. [2] Wang HG, Guan YH, Chen TL, Zhang JT. A study of thermal stresses during laser quenching. J Mater Proc Technol 1997;63:550–3. [3] Li K, Sheng P. Plane stress model for fracture of ceramics during laser cutting. Int J Mach Tools Manuf 1995;35:1493–506. [4] Elperin T, Rudin G. Thermo-elasticity problem for a multilayer coating-substrate assembly irradiated by a laser beam. Int Commun Heat Mass Transfer 1996;23:133–42. [5] Modest MF. Transient elastic and viscoelastic thermal stresses during laser drilling of ceramics. ASME J Heat Transfer 1998;120:892–8. [6] Schulz D, Becker D, Franke J, Kemmerling R, Herziger G. Heat conduction losses in laser cutting of metals. J Phys D 1993;26:1357–63. [7] Shuja SZ, Yilbas BS. Pulsative heating of surfaces. Int J Heat Mass Transfer 1998;41:3899–918. [8] Neto D, Lima C. Non-linear three-dimensional temperature profiles in pulsed laser heated solids. J Phys D 1994;27:1795–804. [9] Yilbas BS. Laser heating process and experimental validation. Int J Heat Mass Transfer 1997;40(5):1131–43. [10] Mason WP, Thurston RN. Physical acoustics, vol. 18. San Diego: Academic Press, 1988. p. 21–123. [11] Scruby CB, Drain LE. Laser ultrasonics; techniques and applications. Bristol, UK: Adam Hilger, 1990. [12] Hsu NN. US National Bureau of Standarts Report, NBSIR85-3234, 1985. [13] Auld BA. Acoustic fields and waves in solids, vol. 1. Florida: Krieger, 1990. [14] Dubois M, Enguehard F, Bertland L. Modeling of laser thermoelastic generation of ultrasound in an orthotropic medium. Appl Phys Lett 1994;65:554–6. [15] Cheng J, Wu L, Zong S. Thermoelastic response of pulsed photothermal deformation of thin plates. J Appl Phys 1994;76:716–22. [16] Graff KF. Wave motion in elastic solids. New York: Oxford University Press, 1975. [17] Dewhurst RJ, Edwards C, Mckie ADW, Palmer SB. Estimation of thin metal sheet using laser generated ultrasound. Appl Phys Lett 1987;51(14):1066–8. [18] Cheng J, Bertelot Y. Theory of laser generated transient Lamb waves in orthotropic paltes. J Phys D 1996;29:1857–67. [19] Shuja SZ, Yilbas BS. Gas-assisted repetitive pulsed heating of a steel surface. Proc Inst Mech Eng, Part C 1998;212:741–57. [20] Gotski TB, Hussaini MY, Lumley JL. Simulation and modeling of turbulent flows. New York: Oxford University Press, 1996.
B.S. Yilbas et al. / Optics and Lasers in Engineering 37 (2002) 63–83
83
[21] Incropera FP, Dewitt DP. Introduction to heat transfer F appendix. New York: Wiley, 1985. p. 667–96. [22] Versteeg HK, Malalasekera W. An introduction to computational fluid dynamics: the finite volume method. New York: Longman Scientific and Technical, 1995. [23] Shuja SZ. Simulation of three-dimensional gas-assisted heating of solid substance. Ph.D. thesis, Mechanical Engineering Department, KFUPM, Dhahran, Saudi Arabia. [24] Patankar SV. Computer analysis of fluid flow and heat transfer. Swansea: Pinridge, 1981. p. 223–52 [Chapter 8]. [25] Bathe KJ. Finite element procedures. Englewood Cliffs, NJ: Prentice-Hall, 1996. p. 148–92 [Chapter 6]. [26] Doyle JF. Wave propagation in structures. New York: Springer, 1989. [27] Mahmood F. Transient response of varying thickness Timonshenko beam under poind load. M.Sc. thesis, Mechanical Engineering Department, KFUPM, Dhahran, Saudi Arabia.