Estimation of heat conduction losses in laser cutting

International Journal of Machine Tools & Manufacture 39 (1999) 431–458

Estimation of heat conduction losses in laser cutting Joseph M. Prusa*, Girish Venkitachalam, Palaniappa A. Molian Mechanical Engineering Department, Iowa State University, Ames, IA 50011, USA Received 1 July 1996; in final form 17 March 1998

Abstract Numerical models of laser cutting are essential for an improved understanding of the process. In order for the models to closely represent real physics the heat lost by conduction during cutting has to be incorporated into them. This paper outlines the details of a mathematical model that is used to estimate heat conduction losses in laser cutting by employing integral methods to solve the three-dimensional, heatconduction equation. Simple, yet accurate, correlations are presented for the conduction heat loss rate, as well as for characteristic thicknesses of the heat affected zone (HAZ). These variables are functions only of Peclet number (Pe), which may be thought of as a dimensionless cutting speed. The correlations are then used to predict (i) the cutting speed in laser cutting, and (ii) the temperature field in the HAZ. The predictions show an excellent match with experimental results. Cutting speed is predicted using an energy balance at the cutting front. An implicit, nonlinear equation for Pe as a function of Stefan number (Ste) and dimensionless heat of combustion results. Ste gives the ratio of sensible heat to latent heat in the workpiece. Detailed analysis of beam absorptivity, coupled with the dimensionless cutting speed result suggests that the dimensional cutting speed, uo, may be correlated as a unique function of the ratio of incoming laser power to workpiece thickness, Qlas/d. This is corroborated by measured cutting rates over a wide range of sample thicknesses and two beam power settings.  1998 Published by Elsevier Science Ltd. All rights reserved. Keywords: Laser; Cutting; Materials; Processing; Heat transfer; Conduction

Nomenclature b c d

Distance from erosion front to beam center, m Specific heat, 660 J/kg K Workpiece thickness, m

* Corresponding author. 0890-6955/98/$—see front matter  1998 Published by Elsevier Science Ltd. All rights reserved. PII: S 0 8 9 0 - 6 9 5 5 ( 9 8 ) 0 0 0 4 1 - 8

432

J.M. Prusa et al. / International Journal of Machine Tools & Manufacture 39 (1999) 431–458

h hf k m ˙ n* uo wo (x¯,y¯,z¯) A F I Pe Qlas Qoxid Qcond Q*c ond (r¯,␪¯ ,z¯) R (R,␪,Z) Ste T T¯ Tˆ (X,Y,Z) ␣ ␥ ␦¯ ␦ ⑀ ␭ ␩ ␮ ␳ ⌬ ⌬H ⌬T ⌫

Heat transfer coefficient, W/m2 K Latent heat of fusion, 275 kJ/kg Thermal conductivity, 38 W/m K Ablation rate, kg/s Complex refractive index Cutting speed, m/s Beam radius, m Cartesian coordinates (m, m, m) Local absorptivity Dimensionless kerf thickness, ⑀/d Local beam intensity, W/m2 Peclet number, ⑀uo/␣ Laser beam power, W Rate of heat input from oxidation, W Rate of conduction heat loss, W Dimensionless conduction heat loss rate Polar coordinates (m, rad, m) Electric field reflectance of electromagnetic wave Dimensionless polar coordinates Stefan number, c⌬T/hf Dimensionless temperature Temperature of the workpiece, K (or °C) Temperature after vertical integration, K (or °C) Dimensionless Cartesian coordinates Thermal diffusivity, k/␳c (m2/s) Cutting front inclination, rad Boundary layer thickness, m Dimensionless boundary layer thickness Half kerf width, m Wavelength, 10.6 ␮m Local angle of incidence at erosion front, rad Stochiometric ratio of FeO:Fe in the melt Density, 7860 kg/m3 Dimensionless temperature perturbation layer thickness Heat of combustion, 258 MJ/kg mol Sensible temperature rise, (Tm − T⬁) Dimensionless heat of combustion, ␮⌬H/(cMFeO⌬T)

Subscripts b cir m par per

Bottom of cut zone Circularly polarized wave Melting point (Tm = 1810 K) Plane parallel polarized wave (P-polarization) Plane perpendicular polarized wave (S-polarization)

J.M. Prusa et al. / International Journal of Machine Tools & Manufacture 39 (1999) 431–458

t 0 1 2 ⬁

433

Top of cut zone Corresponds to x¯ = 0 = X Region I Region II Reference condition of workpiece

Superscripts ⬘ ⬙

First derivative Second derivative

1. Introduction The role played by lasers in manufacturing is a very significant one. Today lasers are employed for cutting, welding, heat treating, vapour deposition, and a number of other applications. Among all the possible applications of lasers, cutting is the most common one, and it accounts for around two-thirds of the total applications. According to published statistics, the laser cutting industry accounts for nearly US$600 million per annum with a projected annual growth rate of about 25%. Over the years, laser cutting has developed from a very crude form to a highly sophisticated process capable of being interfaced with CAD/CAM tools. However, controlling the laser cutting process by numerical means is not an easy task because the process is primarily of a thermal nature. The behaviour of the process near a corner will then be different from its behaviour in the mid-region of a section. This means that adjustments have to be made in a number of situations to maintain the uniformity of the cut. This calls for the knowledge of what parameters to change, and by how much. In order to generate CNC part programs, important information that is required would include cutting speed, laser power, and gas pressure for a given material and sheet thickness. It is in these respects that a laser cutting model becomes highly relevant. In order for a laser cutting model to predict the cutting speed or other processing parameters, the easiest and most direct method would be to solve the overall energy balance equation. Assuming that no evaporative latent heating occurs, this is: Qlas + Qoxid = m ˙ (c⌬T + hf) + Qcond

(1)

where m ˙ = 2␳⑀duo. Energy is lost from the cutting zone by conduction, convection, and radiation. It has been reported by Lim [1] that the major source of heat loss in laser cutting is conduction, while the contributions of convection and radiation are so small that they can be neglected without any loss of accuracy. This observation has also been supported by Powell [2] and Vicanek and Simon [3]. Thus, for the model to work, the heat loss by conduction as a function of the processing parameters has to be computed explicitly. 2. Literature review Laser cutting models can broadly be classified into two categories—those that set out to predict the processing parameters for a set of input conditions, and those that aim to model the behaviour

434

J.M. Prusa et al. / International Journal of Machine Tools & Manufacture 39 (1999) 431–458

of the process to provide a better understanding of the process. These models focus on a number of aspects of laser cutting which include heat transfer phenomena, fluid dynamics of the melt and gas jet, surface quality, and oxidation reactions. A literature review of modeling of heat transfer phenomena reveals that many studies have been carried out covering both the above categories. The models range from extremely simple ones, employing a number of assumptions, to highly detailed ones that endeavour to model the actual situation as closely as possible. The earliest models [4,5] solved the energy balance Eq. (1) assuming that the heat losses were negligible. It is not surprising to find that such models tend to over-predict the results, especially when used with thicker materials. Another approach adopted by a number of researchers [6,7] was based on the solution of the moving line source heat equation given by Rosenthal [8]. Gonsalves and Duley [9] employed a point source model to predict various processing parameters in cutting thin sheets with a CW laser. Bunting and Cornfield [10] developed a relationship between incident power density per unit thickness and the cut speed as a function of the thermal properties of the material by solving the moving line source equation and considering the laser beam as a cylindrical heat source. Hsu [11] modified Bunting and Cornfield’s [10] model to incorporate the effects of varying power density of the laser beam. Analytical calculations of heat loss by conduction were performed by Schulz et al. [12] by employing a cylindrical heat source model. However, the results are presented as a function of the width of the melting front, a parameter that is difficult to determine in practice. Although it is possible to carry out calorimetric experiments to determine the conduction heat loss during laser cutting [13], experiments have to be carried out for each material, sheet thickness, and processing condition of interest. It is much simpler to have an accurate analytical result that may be used to estimate this loss. Such an analytical result would also be useful in helping to organize experimental data into a simple correlation. 3. Objective of the present study Only a few of the large number of models that have been published are capable of predicting fairly close values of the processing parameters in laser cutting, especially of thicker materials. Furthermore, the results produced by many models are not easily reinterpreted for different materials or sheet thicknesses, which could be a disadvantage in an industrial setting. The main objective of the present study is to analytically calculate the conduction heat transfer loss in laser cutting and present this result in an elementary form. The resulting dimensionless correlation only requires that the material thermophysical properties, cutting speed, and the kerf width be known. In an industrial setting, these are all quite measurable properties. The conduction loss correlation is incorporated into an energy balance at the erosion front. This results in an implicit expression that is used to estimate optimum cutting speed. Our criterion for optimum cutting speed is empirical in nature, and is based upon the fastest sustainable cutting speed that makes a clean cut (no dross, waviness, burning, rewelding, and so on). Once this cutting speed has been established, slightly higher speeds result in no cut, while slower ones produce wider kerfs or inferior quality cuts. What this implies for the model is no evaporative latent heating or other sources of energy loss. Since our model does not incorporate these types of losses, it should—and does—predict cutting speeds that match actual cutting trials. This study also characterizes the heat affected zone (HAZ) both in terms of its size and the temperature field within. These results are shown to

J.M. Prusa et al. / International Journal of Machine Tools & Manufacture 39 (1999) 431–458

435

compare well with thermocouple measurements. Finally, detailed computations of the integrated absorptivity at the erosion front show that it is an extremely sensitive function of the position of the laser relative to the erosion front. This sensitivity provides a mechanism whereby the quality of a cut can change enormously based upon edge effects or the hitting of a surface defect on the workpiece. 4. Model and method of solution Fig. 1 shows the phenomena taking place in the cutting zone during laser cutting and the numerical parameters that are of interest in the present study. The laser beam moves in the x¯ direction generating a kerf having a width 2⑀. The y¯ axis corresponds to the direction perpendicular to the motion of the laser beam along the surface, while z¯ is oriented parallel to the laser beam axis. The material that is melted (by heat supplied by the beam and the oxygen jet) is blown from the cutting zone by the coaxial oxygen jet. Heat lost by conduction from the cutting zone is confined to a thermal boundary layer of thickness ␦¯ which is a function of time and spatial location along the cut. This thermal boundary layer, or penetration depth, represents the heat affected zone in cutting. In order to model the heat conduction effects, the cutting zone is divided into two regions. Region I, the region immediately preceding the beam, is modeled by using the polar coordinate system. Region II, the post-beam region, is modeled using the Cartesian coordinate system, and the initial condition for solving the heat conduction equation in this region is supplied by the solution of Region I. Even though considerable effort has been made to model the real situation in laser cutting, a few simplifying assumptions are still employed. It is assumed that the kerf has constant width. In laser cutting, a variation in kerf width from top to bottom is common. However, cutting experi-

Fig. 1.

Cutting zone in laser cutting.

436

J.M. Prusa et al. / International Journal of Machine Tools & Manufacture 39 (1999) 431–458

ments on various sheet thicknesses using a 1.5 kW CW CO2 laser have shown that this kerf variation is very small under optimum cutting conditions; optimum condition being the set of processing parameters that produces cuts of acceptable quality at the maximum possible speed. This being the case, the error produced by a uniform kerf assumption should be marginal. It is also assumed that the thermophysical properties are independent of temperature. 4.1. Dimensional formulation The heat conduction equation transformed to a moving laser coordinate system can be expressed in Regions I and II as:

再

u⬁ − cos␪

u⬁

冉 冊 冎 再 冉 冊

∂T1 1 ∂2T1 ∂2T1 ∂T¯1 sin␪ ∂T1 1 ∂ = ␣ r + + 2 + ¯ ∂r¯ r¯ ∂␪¯ r¯ ∂r¯ ∂r¯ r¯2 ∂␪¯ 2 ∂z¯

冉

∂T2 ∂2T2 ∂2T2 ∂2T2 + 2 + 2 =␣ ∂x¯ ∂x¯2 ∂y¯ ∂z¯

冎

冊

(2a)

(2b)

The boundary conditions in the two regions are: T1 = Tm at r¯ = ⑀; T1→T⬁ as r¯→⬁; ∂T1/∂␪¯ = 0 at ␪¯ = 0; − (∂T1/∂␪¯ )/r¯ = ∂T2/∂x¯ at (␪¯ ,x¯) = ( ± ␲/2,0); − k(∂T1/∂z¯) = ht(T1 − Tt) at z¯ = d; and − k(∂T1/∂z¯) = hb(T1 − Tb) at z¯ = 0; (Region I) T2 = T1 at x¯ = 0; T2→T⬁ as y¯→ ± ⬁; ∂T2/∂y¯ = 0 at y¯ = ± ⑀; − k(∂T2/∂z¯) = ht(T2 − Tt) at z¯ = d; and − k(∂T2/∂z¯) = hb(T2 − Tb) at z¯ = 0; (Region II) The fourth boundary condition (flux match) for Region I and first boundary condition (temperature match) for Region II serve to couple the solutions in both the regions, while the last two conditions in both the regions consider the convective heat transfer at the top and bottom of the material. The above equations are integrated vertically through the workpiece thickness so that a new two-dimensional temperature, Tˆ, can be defined. This procedure reduces the 3D equations into 2D equations, while still accounting for 3D and surface effects. No approximation is introduced by this averaging process—the results still represent an exact analytical formulation. The vertically averaged temperature profiles and the resulting equations are:

冕 d

冕 d

1 1 T1(r¯,␪¯ ,z¯)dz¯; and Tˆ2(x¯,y¯) = T2(x¯,y¯,z¯)dz¯ Tˆ1(r¯,␪¯ ) = d d 0

0

(3a)

J.M. Prusa et al. / International Journal of Machine Tools & Manufacture 39 (1999) 431–458

再

u⬁ − cos␪

u⬁

冉 冊 冎 再 冉 冊

∂Tˆ1 1 ∂2Tˆ1 ∂Tˆ1 sin␪ ∂Tˆ1 1 ∂ + = ␣ r + + ␤1Tˆ1 + ␤2 ¯ ∂r¯ r¯ ∂␪¯ r¯ ∂r¯ ∂r¯ r¯2 ∂␪¯ 2

冉

∂Tˆ2 ∂2Tˆ2 ∂2Tˆ2 + 2 + ␤1Tˆ2 + ␤2 =␣ ∂x¯ ∂x¯2 ∂y¯

冊

冎

437

(3b)

(3c)

where

␤1 = (hb − ht)/kd and ␤2 = (hbTb − htTt)/kd. The terms ␤1 and ␤2 result from the convective boundary conditions at the top and bottom of the material and thus account for surface convective effects. The heat loss by conduction can be calculated from the relation:

冕

␲/2

Qcond = − 2⑀kd (∂Tˆ1/∂r¯)⑀d␪¯

(4)

0

More generally, the integrand in Eq. (4) should be replaced by ⵜTˆ·dA where dA is the infinitesimal area element perpendicular to the unit normal along non-circular erosion fronts. This generalization is important for highly angled erosion fronts or non-symmetrical intensity distributions, such as the line beam [14]. 4.2. Dimensionless formulation The above formulation is nondimensionalized to: 쐌 generalize the formulation and make it independent of any particular unit system 쐌 frame the equations in terms of dimensionless groups that characterize the physics of the problem 쐌 reduce the number of independent parameters Cutting speed, kerf width and diffusivity are grouped together to give the Peclet number, Pe, which represents the ratio of convection to thermal diffusion speeds. Chain rules are applied to Eq. (3b) and (3c) to yield the dimensionless conduction equations. The dimensionless variables, parameters, as well as equations are:

冉

冊

Tˆ1 − T⬁ ; ␦1 = ␦1/⑀; (Region I) R = r¯/␦¯ 1; ␪ = ␪¯ ; and T1 = Tm − T⬁

438

J.M. Prusa et al. / International Journal of Machine Tools & Manufacture 39 (1999) 431–458

再

Pe − (␦1sin␪)⬘

冉 冊 冎

冉 冊 再

∂T1 1 ∂ ∂T1 1 ∂2T1 ␦1sin␪ ∂T1 = R + 2 − N1(T1) + ∂R R ∂␪ R ∂R ∂R R ∂␪2

+ F2(␥1T1 + ␥2)␦ 21

(5)

where Pe = ⑀u0/␣; F = ⑀/d; ␥1 = ␤1d 2; ␥2 = and N1 =

冉 冊

冉

冎

冉

冊

␤1T⬁ + ␤2 2 d ; Tm − T⬁

冊

冉 冊

冉

冊

∂ 2R␦1⬘ ∂2 ␦1⬙ 2␦⬘21 R␦1⬘ + − 2 R − ␦1 ∂R∂␪ ␦1 ␦1 ∂R ␦1

2

∂2 . ∂R2

Tˆ2 − T⬁ ; ␦2 = ␦2/⑀; (Region II) Y = (y¯ − ⑀)/␦¯ 2; X = x¯/⑀; and T2 = Tm − T⬁

冉

Pe ␦ 22 where N2 =

冊

再

冎

∂T2 ∂2T2 ∂T2 ∂2T2 2 + ␦2␦2⬘ = + ␦ − N2(T2) + F2(␥1T2 + ␥2)␦ 22 2 ∂X ∂Y ∂Y2 ∂X2

冉 冊

冉

冊

冉 冊

∂ 2Y␦2⬘ ∂2 ␦2⬙ 2␦⬘22 Y␦2⬘ + − − 2 Y ␦2 ∂X∂Y ␦2 ␦2 ∂Y ␦2

2

(6)

∂2 . ∂Y2

The ranges of the dependent and independent variables in the above equations are 0 ⱕ T1,T2 ⱕ 1; 1/␦1 ⱕ R ⱕ 1 corresponding ⑀ ⱕ r¯ ⱕ ␦1; -␲/2 ⱕ ␪ ⱕ ␲/2; and 0 ⱕ Y ⱕ 1 corresponding to ⑀ ⱕ y¯ ⱕ ⑀ + ␦2. ␦1, ␦2, and X may all become indefinitely large. The dimensionless Y is correct only for the upper half plane, y¯ ⱖ 0. An evaluation of the terms ␥1 and ␥2 in Eqs. (5) and (6) for typical materials (mild steel) and cutting thicknesses (6.35 mm) results in values of the order of 0.1. In laser cutting, since the kerf width is small compared to the sheet thickness, the term F2 is also very small, typically of the order of 10−3. Thus the last term in Eqs. (5) and (6) can be neglected for all practical cases. Scale analysis reveals that the axial (circumferential) diffusion contribution is much smaller than the transverse contribution for large Pe case (high speed cutting). However, the circumferential diffusion terms in Region I are retained as present day cutting speeds correspond to Pe of order one. The heat equation is strongly elliptic in this region. 4.3. Integral method The mathematical model formulated above in Regions I and II in terms of partial derivatives is solved using an integral method. This method is a semi-analytical technique capable of producing

J.M. Prusa et al. / International Journal of Machine Tools & Manufacture 39 (1999) 431–458

439

approximate solutions to transient heat conduction problems. Unlike other methods, the integral method does not linearize the problem; it reduces the order of the nonlinear boundary value problem, often to an ordinary initial value problem that can be solved relatively easily. The reader is directed to Goodman [15] and Arpaci [16] for more information on integral methods. A review of integral methods applied to melting problems may be found in Yao and Prusa [17]. The two regions will be treated separately for presentation purposes. 4.3.1. Region I In order to reduce the partial differential Eq. (5) to an ordinary differential equation, a quadratic polynomial approximation in radius is used for determining the temperature distribution: T1 = a0 + a1R + a2R2

(7)

This profile can be represented graphically, as shown in Fig. 2. Here a0, a1, and a2 are unknown coefficients and are determined by forcing Eq. (7) to satisfy the dimensionless conditions: T1(1/␦1) = 1, T1(1) = 0, and (∂T1/∂R)1 = 0 It has to be noted that the last two boundary conditions shown above are approximate and follow from the concept of the thermal boundary layer. Thus they replace the condition T1→0 as R→⬁. The temperature profile so obtained is: T1 = ␦ 21

冉 冊 1−R 1 − ␦1

2

(8)

The penetration depth, ␦1(␪) can now be determined by substituting Eq. (8) into Eq. (5)—with the negligible terms of O(F2) (e.g., order of magnitude of F2) being dropped—and then integrating

Fig. 2. Dimensionless temperature profile in Region I.

440

J.M. Prusa et al. / International Journal of Machine Tools & Manufacture 39 (1999) 431–458

(with volume weighting) the resulting equation in the R direction from R = 1/␦1 to R = 1. This procedure eliminates the explicit R dependence and results in the ordinary differential equation connecting ␦1 and ␪ as shown below: 兵3␦1(1 − ␦1)(1 − ␦ 21 + 2␦1ln␦1)其␦1⬙ + 兵3␦1␦1⬘[(1 − ␦1)(5 + ␦1) + 2(1 + 2␦1)ln␦1] − Pe␦1(1 − ␦1)4sin␪其␦1⬘ − 兵(1 − ␦1)3[6 + Pe(1 − ␦1)(2 + ␦1)cos␪]其␦1 = 0

(9)

Primes denote differentiation with respect to ␪. The above equation can be solved for a particular Peclet number to obtain the extent of the heat affected zone corresponding to that particular set of processing parameters. The solution procedure for the full equation will be discussed in the next section. A dimensionless heat loss Q*cond can be defined that is connected to the dimensional heat loss as follows: Q*cond = Qcond/(kd⌬T)

(10)

The dimensionless heat loss (for semi-circular erosion fronts) can be calculated from the relation

冕

␲/2

冕

␲/2

Q*cond = − 2 ␦ 1− 1(∂T1/∂R)1/␦1d␪ = 4 (␦1 − 1)−1d␪ 0

(11)

0

The heat loss by conduction can thus be calculated by employing Eqs. (10) and (11) against different Peclet numbers which represent the grouping of the process parameters. The primary advantage of the above procedure is that, once the dimensionless heat transfer is calculated, conduction loss in any material and for any thickness can be calculated by employing Eq. (10). Zero order (constant) asymptotic solutions for Eqs. (9) and (11) valid in the limit Pe→⬁ are (␦1 − 1)→2/Pe and Q*cond→␲Pe. First order (showing ␪ dependency), parabolic solutions will be examined later in Section 5. 4.3.2. Region II Even though the basic approach to the solution in Region II is the same as in Region I, there are some technical differences in how the integral method is implemented. A little consideration will reveal that the solution in Region II will have to be matched with that of Region I, and this matching has to be done at the boundary of these two regions. This matching condition is in fact the fundamental reason for considering Region II at all: it provides a boundary condition for the elliptic problem in Region I. The two regions are coupled by noting that the penetration depth and temperature distribution will be the same in both regions at the matching boundary. From Fig. 1 it follows that:

␦1(␪¯ = ␲/2) = ␦2(x¯ = 0) + ⑀ Transforming this equation into dimensionless form:

(12)

J.M. Prusa et al. / International Journal of Machine Tools & Manufacture 39 (1999) 431–458

␦1(␪ = ␲/2) = ␦2(X = 0) + 1

441

(13)

The temperature profile at the Region I/II interface, ␪ = ␲/2 or X = 0, can be obtained from Eq. (8) by substituting the value of ␦1(␪ = ␲/2) into it. The resulting equation is transformed wholly into Cartesian coordinates by employing the relations: x¯ = − r¯cos␪¯ , y¯ = r¯sin␪¯ , and r¯2 = x¯2 + y¯2 A conversion into dimensionless form yields the necessary relation: R2 = 兵X2 + (1 + ␦2Y)2其/␦ 21

(14)

Substituting Eqs. (13) and (14) into Eq. (8) and noting that X = 0 at the common boundary between Regions I and II yields: T2,0 = (1 − Y)2

(15)

which provides the initial condition for the Region II solution. Note that as X increases from its initial value of zero, this initial temperature profile rapidly falls out of step with the new adiabatic boundary condition at the kerf, (∂T2/∂Y)0 = 0. Following the disturbance method of Goodman [15], a correction or perturbation temperature is added to the above profile to satisfy this requirement. This perturbation is added to Eq. (15) to obtain a total temperature profile that will satisfy the boundary conditions. Associated with this perturbation, a second penetration depth ⌬ is defined with the initial condition ⌬ = 0 at X = 0 (⌬ approaches a maximum value of 1 as x¯,X→⬁). Defining the profile in Eq. (15) as the outer profile and the perturbation that has to be added as the inner profile, the total profile can be determined by noting the requirement: ∂T2 ∂Y

|

⬅0=

0

∂T2,out ∂Y

|

0

+

∂T2,in ∂Y

|

(16)

0

A quadratic polynomial distribution, shown in Fig. 3, is also assumed for the inner profile. The associated unknown coefficients are then determined by applying the condition shown in Eq. (16) and the following two boundary conditions: ∂T2,in/∂Y = 0, and T2,in = 0 at Y = ⌬ (corresponding to y¯ = ⑀ + ⌬␦2). The total temperature distribution is then obtained as: T2 =

冦

(1 − Y)2 − ⌬(1 − Y/⌬)2 for 0 ⱕ Y ⱕ ⌬ (1 − Y)2 for ⌬ ⱕ Y ⱕ 1 0 for 1 ⱕ Y

(17)

442

J.M. Prusa et al. / International Journal of Machine Tools & Manufacture 39 (1999) 431–458

Fig. 3.

Dimensionless temperature profile in Region Il.

This total temperature profile is then substituted into a parabolized version (e.g., the diffusion terms in the X direction are dropped) of the dimensionless conduction Eq. (6), in which terms of O(F2) are also dropped. The result is integrated in the Y direction throughout the heat affected zone, from 0 ⱕ Y ⱕ 1. Next, the perturbation temperature, − ⌬(1 − Y/⌬)2, itself is substituted into the same PDE, which is then integrated only in the inner region, 0 ⱕ Y ⱕ ⌬. The following two differential equations are obtained for ␦2, and ⌬: (2␦ 22⌬⌬⬘ + ␦2␦2⬘⌬2)Pe = 6 (perturbation temperature)

(18)

2␦2⌬⌬⬘ + ␦2⬘(⌬2 − 1) = 0 (full energy balance)

(19)

Eqs. (18) and (19) can be solved simultaneously to obtain analytical solutions for the boundary layer profiles, ␦2 and ⌬. The solutions are:

␦2 = (␦ 290 + 12X/Pe)1/2

(20)

⌬ = (1 − ␦90/␦2)1/2

(21)

Here ␦90 is the initial value of ␦2 at X = 0. It is given by the solution from Region I applied to the matching boundary at ␪ = ␲/2 (e.g., at 90°). 5. Numerical method and computational procedure Inasmuch as Region II was modeled primarily to provide a boundary condition for Region I, it sufficed to approximate Region II with a parabolic model. The resulting differential equations

J.M. Prusa et al. / International Journal of Machine Tools & Manufacture 39 (1999) 431–458

443

from the integral method, although nonlinear, were analytically tractable. Unfortunately, the elliptic equation resulting in Region I (Eq. (9)) is not so amenable to an exact solution so a numerical method was employed for solving it. We used an implicit, second order Gauss–Seidel method. In order to obtain an initial profile for the elliptic iterations, a parabolic formulation (in which the circumferential diffusion terms B1 and B2 below were dropped) of the ODE was solved using a second order Runge–Kutta technique. Eq. (9) can be represented as: B1␦1⬙ + (B2 + B3)␦1⬘ + B4␦1 = 0

(22)

where B1 = 3␦1(1 − ␦1)(1 − ␦ 21 + 2␦1ln␦1), B2 = 3␦1␦1⬘[(1 − ␦1)(5 + ␦1) + 2(1 + 2␦1)ln␦1], B3 = − Pe␦1(1 − ␦1)4sin␪, and B4 = (1 − ␦1)3[6 + Pe(1 − ␦1)(2 + ␦1)cos␪] By applying central difference relations to the above equation and rearranging the terms, the following finite difference formula for solving ␦1 was obtained:

冉

冊

冉

冊

冉

冊

2B1 B1 B2 + B3 k B1 B2 + B3 k + 1 + B4 ␦ k1,i+ 1 = + ␦ 1,i + 1 + − ␦ 1,i − 1 2 2 ⌬␪ ⌬␪ 2⌬␪ ⌬␪2 2⌬␪

(23)

where i represents the spatial location of the node, while k stands for the iteration level. The penetration depth, ␦1, in the terms B1 through B4 was evaluated as ␦ k1,i. Since Eq. (9) is second order in ␪, two boundary conditions are required to solve it. One boundary condition comes from the symmetry condition at the leading edge of the boundary layer, (∂T1/∂␪)0 = 0. This condition was translated into finite difference form as:

␦ k1,1+ 1 = (4␦ k1,2 − ␦ k1,3)/3

(24)

The second boundary condition was derived from the matching requirement between Regions I and II. Since the heat flux through the interface from one region to the other must be the same, the flux relations were equated to obtain the necessary condition. Furthermore, as an analytical solution exists in Region II for the penetration depth, this boundary condition was obtained readily. The heat flux relations are: ␦¯ 90

⫺

冕

␦¯ 90

r−1(∂Tˆ1/∂␪¯ )␲/2dr¯ ⫽

⑀

冕

(∂Tˆ2/∂x¯)0dy¯

(25)

⑀

where ␦¯ 90 is the dimensional penetration depth at the matching boundary between Regions I and II. Eq. (25) was transformed into dimensionless form and then integrated incorporating the parabolic temperature profiles determined previously. Finally, Eq. (20) was applied to eliminate explicit dependence on ␦2 and its derivative, yielding the second boundary condition for Region I:

␦ 1⬘ = −

冉 冊再

冋

2 2␦ 21 (␦1 − 1)(1 − 3␦1) ␦1 1+ ln␦1 + 3 Pe ␦1 − 1 (1 − ␦1) 2␦ 21

册冎

−1

(26) ␲/2

444

J.M. Prusa et al. / International Journal of Machine Tools & Manufacture 39 (1999) 431–458

This boundary condition was finite differenced using the same, one-sided, second order accurate form as was used for the symmetry condition, Eq. (24). All ␦1 terms on the RHS of Eq. (26) were evaluated using their previously iterated values. The coding for the finite difference problem was performed using FORTRAN 77. A convergence criterion of |␦ k1,i+ 1 − ␦ k1,i|/␦ k1,i ⱕ 10−6 was employed for the iterations. A local cubic curve smoothing was used on the first few nodes near the symmetry point in the parabolic initial condition for Region I. This smoothing profile simultaneously satisfied the symmetry boundary condition ␦1⬘(0) = 0, a further boundary condition for ␦1⬙(0) given by Eq. (9) applied at ␪ = 0, and an averaged applied fit (in a least squares sense) to the ODE in the small initial interval in which the smoothing was used. This was found to be helpful in reducing the number of elliptic iterations as the parabolic problem is singular at ␪ = 0. The angular resolution used was varied according to the value of Pe so that the Q*cond and ␦90 values computed were converged to at least 3 significant figures. The requisite number of nodes varied from 21 at Pe = 0.1 to 81 at Pe = 10.

6. Results and discussion A series of numerical computations of Eq. (23) were carried out over the Peclet number range, 0.1 ⱕ Pe ⱕ 10.0. This range encompasses the lower Pe values typical of modern laser cutting processes, as well as a higher range that hopefully will be appropriate for future, improved cutting processes. We begin by summarizing the basic predictions of the model for the size of the heat affected zone (HAZ), the conduction heat transfer loss into the workpiece, and some elementary, dimensionless correlations for these predictions. Next, the model results are incorporated into an energy balance for the erosion front in order to predict cutting speeds. A detailed analysis is made of the local absorptivity at the erosion front in order to estimate the average absorptivity during cutting. Finally, comparisons are made with experimental measurements of cutting speeds and temperature fields. 6.1. The HAZ Figs. 4–6 show the dimensionless penetration depths obtained from Eqs. (20), (21) and (23). The results are shown for a Pe range of 0.3 ⱕ Pe ⱕ 10.0. The outer penetration depths, ␦1, and ␦2 decrease with an increase in Pe because this is equivalent to an increase in cutting speed. These penetration depths are summarized in Table 1, which gives the dimensionless HAZ thickness, (␦1 − 1), at the leading edge and at ␪ = 90°. The dimensionless inner, or perturbation penetration depth, ⌬, obtained from Eq. (21) is shown in Fig. 6. ⌬ rapidly increases from its initial value of zero at X = 0 to a maximum value of unity. This means that the change in boundary condition along the cut surface (from the melting temperature to an adiabatic condition) quickly propagates out into and dominates the HAZ in Region II. These results may immediately be translated into physical dimensions by multiplying by the kerf half-thickness, ⑀, according to the dimensionless variable definitions accompanying Eqs. (5) and (6). It is important to note that ␦1 and ␦2 are defined slightly different (recall Eq. (13)). In particular:

J.M. Prusa et al. / International Journal of Machine Tools & Manufacture 39 (1999) 431–458

445

Fig. 4. Dimensionless thermal penetration depth, ␦1, in Region I as a function of ␪ and Pe(0 ⬍ ␪ ⱕ ␲/2 radians).

Fig. 5. Dimensionless thermal penetration depth, ␦2 in Region II as a function of X and Pe.

r¯/⑀ = R␦1 and Y* ⬅ y¯/⑀ = 1 + Y␦2 For example, consider the cutting of a 6.35 mm (0.25⬙) thick mild steel plate at 1000 W. Experiments yielded a cutting speed of 18.6 mm/s (44 ipm) and a kerf width of 0.50 mm (0.0197⬙). Assuming k = 38 W/m K and c = 660 J/kg K—typical integral averaged values (from 300 to 1800 K) for mild steel—this translates into a Pe of 0.64. The dimensional penetration depths

446

J.M. Prusa et al. / International Journal of Machine Tools & Manufacture 39 (1999) 431–458

Fig. 6. Dimensionless perturbation depth, ⌬, in Region II as a function of X and Pe (asymptotic value is 1 for all Pe). Table 1 Numerical solutions for the dimensionless thermal boundary layer thickness Pe

␦1(0) − 1

0.1 0.2 0.3 0.6 0.8 1.0 2.0 3.0 4.0 6.0 8.0 10.0

13.0 6.56 4.47 2.38 1.85 1.52 0.828 0.577 0.446 0.306 0.234 0.189

␦90 26.0 13.7 9.58 5.35 4.24 3.55 2.09 1.55 1.26 0.943 0.770 0.659

computed using this Pe are shown in Figs. 7 and 8. The figures show that the penetration depth in laser cutting is indeed very small, ranging from 0.6 mm immediately in front of the laser (␦1(0) − 1 = 2.25), to 1.3 mm at the boundary between Regions I and II (␦90 = 5.07). In Region II, some 25 mm behind the laser, the outer penetration depth has expanded to 11 mm (␦2 = 44). These estimates for penetration depth overstate the region of high temperature, which is significantly smaller. In particular, the region with a temperature above 200°C extends only 0.4 mm in front of the laser, and is 0.8 mm thick at ␪ = 90°. This isotherm then curves back into the kerf and ends just 1.3 mm downstream from the laser (see Fig. 9). This behaviour is due to the perturbation temperature which is acting so as to bring the total temperature profile into compliance with the adiabatic boundary condition along the kerf for X > 0. The extent of the region

J.M. Prusa et al. / International Journal of Machine Tools & Manufacture 39 (1999) 431–458

447

Fig. 7. Dimensional thermal penetration depth, ⑀(␦1 − 1) (outer edge of HAZ), in Region I for case Pe = 0.64 (cutting of 6.35 mm mild steel plate at 18.6 mm/s (44 ipm) with kerf width of 0.50 mm).

Fig. 8. Dimensional thermal penetration depths, ⑀␦2 (outer edge of HAZ) and ⑀⌬␦2 (outer edge of inner perturbed region), in Region II for case Pe = 0.64 (cutting of 6.35 mm mild steel plate at 18.6 mm/s (44 ipm) with kerf width of 0.50 mm).

over which this perturbation acts is also shown in Fig. 8—it rapidly increases from zero to nearly the same thickness as the outer penetration depth. Fig. 9 shows the isothermal contours in the HAZ at 30, 200, 500, and 800°C, for the case Pe = 0.64. The 30°C isotherm is the edge of the HAZ. The kerf is denoted in bold along the bottom of the figure with the laser beam centred at (X,Y*) = (0,0). The dotted grid shows 1.00 mm squares. The isotherm at 800°C has two branches

448

J.M. Prusa et al. / International Journal of Machine Tools & Manufacture 39 (1999) 431–458

Fig. 9. Isotherms in the HAZ for case Pe = 0.64 (cutting of 6.35 mm mild steel plate at 18.6 mm/s (44 ipm) with kerf width of 0.50 mm). The 30°C isotherm is the outer edge of the HAZ. The heavy contour along the bottom is the edge of the kerf; for X ⬍ 0, this is the erosion front and is at 1535°C. The dotted rectangles denote 1.0 mm squares.

in Region II. The solid downward branch is the parabolic solution used for Region II. Most likely, it shows the isotherm curving too sharply into the kerf because the Region II model does not incorporate diffusion in the X direction. The upper, dashed, branch corresponds to the solution in which ⌬¿1, that is, the kerf is close to the melt temperature for all X. This solution certainly gives an extreme upper bound for the hotter isotherms (it never intersects the kerf). The true elliptic solution should lie somewhere between these extremes. For isotherms that lie further from the laser beam, the ellipticity of Region II becomes less important, e.g., the 200°C isotherm should be more correctly placed than the 800°C isotherm. Ellipticity in Region II also becomes less important for larger values of Pe. 6.2. Conduction loss Fig. 10 shows the dimensionless conduction heat transfer loss, Q*cond, obtained from the numerical solution generated by Eq. (23). The parabolic solution used as the initial condition to begin the elliptic iterations of Eq. (23) is also shown in the figure. Both results for Q*cond are also given in Table 2 as well as the percent difference in heat transfer between the two cases. The table clearly shows that for present day industrial cutting speeds, for which typically 0.3 ⱕ Pe ⱕ 1.5, that the parabolic solution is off by 50–10%. This demonstrates the necessity of using an elliptic formulation as in the present model. For the cutting example considered above, Pe = 0.64, the

J.M. Prusa et al. / International Journal of Machine Tools & Manufacture 39 (1999) 431–458

449

Fig. 10. Dimensionless conduction heat loss rate and enthalpy flux (for mild steel) as a function of Pe. The solid line gives the full elliptic solution for Q*cond (see Table 2) which is approximated by Eq. (30). The dashed line gives the parabolic result for Q*cond that occurs if the coefficients B1 and B2 in Eq. (22) are dropped. The dotted line for H*t is given by Eq. (27). Table 2 Numerical solutions for the dimensionless conduction heat loss during laser cutting Pe 0.1 0.2 0.3 0.6 0.8 1.0 2.0 3.0 4.0 6.0 8.0 10.0

Q*cond 0.394 0.770 1.125 2.090 2.680 3.242 5.831 8.240 10.56 15.06 19.46 23.78

Parabolic 0.86 1.30 1.67 2.63 3.21 3.76 6.29 8.69 10.98 15.54 19.81 24.12

% Change + 120 70 50 25 20 16 7.9 5.5 4.0 3.2 1.8 1.4

numerical integration yields Q*cond = 2.16. In dimensional form this becomes 780 W (via Eq. (10)). This conduction loss is about 30% greater than the enthalpy flux, Ht = m ˙ (c⌬T + hf), of 590 W, which is the power needed to (i) sensibly heat the steel to the melting point and then (ii) melt it. Fig. 10 also shows the dimensionless enthalpy flux, given by: H*t ⬅ Ht/(kd⌬T) = 2Pe(1 + 1/Ste)

(27)

where Ste ⬅ c⌬T/hf is the Stefan number, and is a measure of the ratio of sensible to latent

450

J.M. Prusa et al. / International Journal of Machine Tools & Manufacture 39 (1999) 431–458

heating required to melt the workpiece [17]. For mild steel, Ste ⬇ 3.60 so H*t ⬇ 2.56Pe. While the conduction loss and the enthalpy flux remain comparable for all Pe shown in the figure, the conduction loss is larger for smaller Pe, but eventually becomes somewhat smaller for larger Pe. Fig. 11 shows the dimensional conduction and enthalpy losses as a function of the cutting speed for mild steel with a thickness of 6.35 mm, assuming an average kerf width of 0.50 mm. 6.3. Correlations for heat transfer and thermal penetration depth The following correlations are found to duplicate the solutions given in Tables 1 and 2 with a maximum error of ± 3% within the specified Peclet number ranges:

␦0 = 1.53Pe−0.905 for 0.15 ⱕ Pe ⱕ 10

(28)

␦90 = 3.70Pe−0.769 for 0.25 ⱕ Pe ⱕ 10

(29)

Q*cond = 3.20Pe+0.868 for 0.20 ⱕ Pe ⱕ 10

(30)

In correlation (28), ␦0 denotes ␦1(0) − 1. The maximum errors all occur at the lowest value of Peclet number cited. The average error over most of the Pe ranges is closer to ± 1.5%. 6.4. Experimental verification Integral method solutions typically have errors on the order of 1–10% [15]. In order to verify the accuracy of the model, two types of comparisons were made. First, the heat conduction results

Fig. 11. Dimensional conduction heat loss rate and enthalpy flux as a function of cutting speed for 6.35 mm mild steel plate.

J.M. Prusa et al. / International Journal of Machine Tools & Manufacture 39 (1999) 431–458

451

were incorporated into the energy balance Eq. (1) to predict the cutting speed for mild steel of different thicknesses. In the second comparison, the model was used to predict the temperature field in the workpiece in the neighbourhood of the laser beam. These predictions were compared with experimental results obtained with a 1.5 kW CO2 laser. The laser delivers a raw beam of diameter 0.019 m which is focused to a spot of diameter 0.20 mm by a 0.127 m focal length ZnSe lens. The beam was circularly polarized and oxygen was used as the assist gas. 6.4.1. Cutting speed predictions Inspection of the energy balance, Eq. (1), shows that three heat transfer terms, Qlas, Qoxid and Qcond must be estimated if it is to be used to determine the cutting speed. The conduction loss, Qcond is estimated from the results already presented in this study using Eq. (30) and Eq. (10). A caveat, of course, is that Eq. (30) requires the Pe as input and this is none other than the dimensionless cutting speed. Thus Eq. (1) gives an implicit function for the cutting speed. The total power absorbed by the laser beam, Qlas, depends on the local absorptivity, A, and the inclination angle, ␥, of the cutting front. Fig. 12 shows the top view of the cutting kerf detailing some of the parameters used in evaluating the absorbed laser power. o is the offset distance between the kerf centre and the beam centre, and is a relatively unknown parameter. In the present model, calculations were performed using different values of o and it was found that the absorbed power was extremely sensitive to the exact value of o. Typically, a value of o = 0.1 mm resulted in maximum absorptivities on the order of 0.3–0.45. These maxima were extremely sensitive functions of ␥, which was itself found to be a function of workpiece thickness. Increasing the offset distance by as little as 0.1 mm could result in sharp drops down to the normal absorptivity value—the minimum possible (typically 0.2–0.25). This very sharp, non-linear response may easily cause a cut that is proceeding nicely to degenerate into a poor cut (or none at all) if a perturbation to the erosion front occurs such that the effect is a tiny increase in offset distance. Since the model was to simulate optimum cutting conditions, we chose an offset distance of 0.1 mm (maximum absorptivity) for the following predictions. The absorbed beam power was computed from:

Fig. 12. Top view of beam/erosion front geometry.

452

J.M. Prusa et al. / International Journal of Machine Tools & Manufacture 39 (1999) 431–458

Qlas ⫽

冕冕

AcirI(b)cot␩da

(31)

a

Here b is the radial distance of the infinitesimal area element cot␩da along the erosion front from the beam axis. This area element is the projected area perpendicular to the beam axis as ␩ is the local angle of incidence. Assuming that the erosion front is an oblique, semi-circular cylinder at inclination ␥, then da = ⑀dz¯d␾ and ␩ and b may be written as the following functions of ␾ and z¯, along with the intensity profile for a Gaussian beam (TEM00 mode):

␩ = cos−1

冤

sin␾ √sin ␾ + tan2␥ 2

冥

(32a)

b = {⑀2 + [o + (d − z¯)cot␥]2 − 2⑀[o + (d − z¯)cot␥]sin␾}1/2

(32b)

I(b) = I0exp( − 2b2/w20)

(32c)

The local absorptivity is a highly sensitive function of the angle of incidence as well as the type of polarization. In particular, plane parallel polarized waves are completely absorbed at the Brewster angle (about 86.0° in the present study). Thus erosion front inclination angles very near ␥ ⬇ 86.0° (this angle should be characteristic of the fastest cutting speeds) have the greatest overall absorption (for plane parallel and circularly polarized beams). The local absorptivity drops off on either side of the Brewster angle, going to zero at 90° and decaying more slowly to the normal absorptivity at 0°. The local absorptivity also depends upon the refractive index of the material, temperature, and presence of oxide layer. Petring et al. [18] provide relationships of absorptivity of plane polarized laser beams as a function of refractive index and angle of incidence. Since the laser beam in our study is circularly polarized, its electric field reflectance, R, is computed by integrating the vector sum of the plane parallel and perpendicular polarized wave reflectances over all polarization angles [19]. Given that A = 1 − R2, the final result is: Acir = (Apar + Aper)/2, where Apar = 1 −

冋

Aper = 1 −

冋

and

(33a)

n*2cos␩ − (n*2 − sin2␩)1/2 n*2cos␩ + (n*2 − sin2␩)1/2 cos␩ − (n*2 − sin2␩)1/2 cos␩ + (n*2 − sin2␩)1/2

册

册

2

(33b)

2

(33c)

J.M. Prusa et al. / International Journal of Machine Tools & Manufacture 39 (1999) 431–458

453

where, n* = nr + ini is the complex refractive index and i is √ − 1. Here n*2 denotes the square of the modulus of n*. The value of n*2 is a function of numerous melt/assist gas thermophysical and electromagnetic properties. The value computed for the present study was n*2 ⬇ 250 [18,20– 22]. Given this value of refractive index, the maximum local absorptivity (at the Brewster angle) was 0.51, while the normal absorptivity was 0.22. Eq. (31) can be readily integrated after substituting Eq. (32a), (32b) and (32c) and Eq. (33a), (33b) and (33c) into it. The resulting integral was evaluated using the trapezoidal rule over 0 ⱕ ␾ ⱕ ␲/2 and 0 ⱕ z¯ ⱕ d and multiplied by two. A resolution of 31 × 31 area elements was found be sufficient to reduce trucation error to 0.1% or less. The integrated absorptivities were then correlated with workpiece thickness (note that ␥ is an input parameter). Data for this correlation came from direct measurements of ␥ using a scanning electron microscope, as well as from several reversed energy balance computations that were used to establish the integrated absorptivity. Kerf thicknesses were measured using an optical microscope and feeler gauges. It is noteworthy that the range of integrated absorptivities was found to be 0.2–0.45, in good agreement with the electromagnetic field theory prediction. These integrated values are summarized in Table 3. The first case, with d = 0.91 mm, used Ar as the assist gas (no combustion). The beam power was 1.5 kW for the two thickest cases, with d = 9.53 mm, and 1.0 kW for all the others. While the data are rough, they are of sufficient quality that when combined with the field theory, a reasonable correlation can be developed for the integrated absorptivity of a circularly polarized beam: Aint ⬇ 0.015 + 0.245Tan−1(1.25d) for 0.9 ⱕ d ⱕ 12.7 mm

(34)

We note that for given ␥, as d→0, less of the beam is intercepted by the workpiece and so Aint→0 (beam diameter becomes large compared to d). To partially compensate for this geometrical effect, ␥ must decrease as d→0, ensuring that the workpiece intercepts enough of the beam for cutting to occur. The last term in the energy balance, Qoxid, is due to combustion of the iron with oxygen in the gas assist jet. It is calculated from [23]: Qoxid = ␮m ˙ ⌬H/MFeO

(35)

where ␮ is the stochiometric ratio of FeO:Fe in the melt and is assumed equal to 0.5, ⌬H = 258 Table 3 Experimental cutting parameters (TEM00 mode) d (mm)

uo(ms−1)

⑀(mm)

0.89 1.45 1.45 3.17 6.37 6.35 9.53 9.53

0.0762 0.106 0.123 * 0.0169 * * *

0.16 0.13 0.14 0.21 0.31 0.25 0.40 0.42

␥ (degrees) * * * 70 * 80 75 82

Aint 0.25 0.17 0.20 0.31 0.34 0.41 0.35 0.45

454

J.M. Prusa et al. / International Journal of Machine Tools & Manufacture 39 (1999) 431–458

MJ/kg mol is the heat of combustion, and MFeO = 71.8 kg/kg mol is the molecular weight of FeO. Note that Eq. (35) contributes another implicit cutting speed term to the overall energy balance. Eq. (1) may now be completely rewritten in dimensionless form as: AintQ*las = 2Pe(1 + 1/Ste − ⌫) + 3.20Pe0.868 * las

(36a) * las

where AintQ is the dimensionless absorbed laser power, with Q being defined as in Eq. (10). The last term on the RHS of Eq. (36a) is the conduction loss, given by Eq. (30). The enthalpy flux, Eq. (27) may also be recognized in Eq. (36a). The remaining term, − 2Pe⌫, is the dimensionless version of Eq. (35), where ⌫ = ␮⌬H/(c⌬TMFeo) is the ratio of heat of combustion to sensible heat in the melt. For the cutting of mild steel using an oxygen gas assist jet, ⌫ ⬇ 1.80. Thus Eq. (36a) reduces to: AintQ*las = − 1.05Pe + 3.20Pe0.868 (mild steel, O2 jet)

(36b)

Eq. (36a) and (36b) are noteworthy because they predict that for a given material and assist jet gas, that the dimensionless cutting speed, Pe, is solely a function of the dimensionless absorbed laser power, Pe = Pe(AintQ*las). Furthermore, given the definition of Pe following Eq. (5), if the integrated absorptivity and kerf thickness can be correlated as a function of workpiece thickness, then dimensional cutting speed, uo, can be reduced to a function of only one variable—the ratio of laser beam power to workpiece thickness, Qlas/d. Eq. (34) presents the integrated absorptivity correlation. Using the results of Hsu and Molian [24], the following correlation for kerf thickness was developed: ⑀(d) = 0.091 + 0.0127d (with an average error of ± 3% in the range 1.0 ⱕ d ⱕ 12.7 mm, with Qlas set to 1.0 or 1.5 kW). Fig. 13 shows the dimensional relation, uo = uo(Qlas/d), as well as a comparison between the

Fig. 13. Comparison of predicted and experimental cutting speeds for mild steel with oxygen assist gas. The dotted lines show the effects of uncertainty of ± 10% variation in the values of c and k compared to the nominal values. The dashed line shows a case of constant integrated absorptivity (at the nominal c and k values). The experimental range in workpiece thickness, d, is from 1.26 mm at 1.0 kW (on the right) to 12.7 mm at 1.5 kW (on the left).

J.M. Prusa et al. / International Journal of Machine Tools & Manufacture 39 (1999) 431–458

455

predicted cutting speeds and those observed in Hsu and Molian [24]. The data for Qlas/d range from the cutting of 12.7 mm material at 1.5 kW (on the left) to the cutting of 1.26 mm material at 1.0 kW (on the right). The solid line shows our result using the nominal values of k and c, and Eq. (34). The comparison with the experimental results is remarkably good. The effects of unknown variations in thermophysical properties are shown by the two dotted lines above and below the nominal solution. These solutions correspond to ± 10% variations in the values of k and c such that the variations in the solutions constructively add. Finally, the dashed line in Fig. 13 shows the predicted solution when a constant integrated absorptivity of 0.37 is used rather than Eq. (34). Although this solution is good for thicker materials, 0 ⱕ Qlas/d ⱕ 0.3 kW/mm, it gives too high a cutting speed for thinner materials, Qlas/d > 0.3 kW/mm. Clearly, Aint must decrease as the material gets thinner. 6.4.2. Temperature field predictions Fig. 14 shows the temperature time histories (TTHs) recorded by two type K thermocouples during a cutting trial on a 9.35 mm (nominally 3/8⬙) mild steel plate as a set of data points. The thermocouples were spot welded to the surface, centred at 5.9 and 15.5 mm from the kerf edge, and had beads approximately 1 mm in diameter. The cutting speed was 1.69 cm/s (40 ipm) with a beam power of 1550 W and a donut, or TEM01 beam mode. The kerf had an average thickness of 0.439 mm (0.0173⬙), resulting in a Peclet number of Pe = 0.508. The Region II integral solution for temperature at these locations is also shown plotted in the figure, in the form of a solid line bounded above and below by dotted lines. The solid curves correspond to the solution using the nominal property and bead position values listed above. The dotted curves above the nominal solutions correspond to 10% variations in k and c such that the Peclet number is decreased to Pe

Fig. 14. Comparison of predicted and experimental temperature time histories during the cutting of 9.35 mm mild steel plate. The nominal cases refer to predictions based upon the nominal c and k values, and thermocouple junction position at the centre of the thermocouple bead. The dotted lines above and below the nominal cases show effects of uncertainty of ± 10% variation in properties and ± 0.5 mm variation in position of the thermocouple junctions.

456

J.M. Prusa et al. / International Journal of Machine Tools & Manufacture 39 (1999) 431–458

= 0.416, and y¯ − ⑀ = 5.4 and 15.0 mm, or 0.5 mm less than the nominal values. The dotted curves below the nominal solutions correspond to 10% variations in k and c such that the Peclet number is increased to Pe = 0.621, and y¯ − ⑀ = 6.4 and 16.0 mm, or 0.5 mm more than the nominal value. The variations in properties and y¯ − ⑀ are chosen so that they constructively add and thus produce the maximum variation in TTH from the nominal solution. The integral solution shows the correct qualitative behaviour, and has the correct temperature maximum. But the thermocouple readings at y¯ − ⑀ = 5.9 mm show the TTH rising slowly above the initial plate temperature, T⬁ = 24°C, before the integral solution’s meteoric rise and then reaching a maximum temperature somewhat later than the integral solution. The readings at y¯ − ⑀ = 15.5 mm show qualitatively similar behaviour. These discrepancies are most likely due to the parabolic nature of the Region II solutions; the full elliptic nature of the Region I solution manifests itself in the TTHs only as an initial condition for ␦2 as given by correlation (29). 7. Conclusions An integral method was used to develop a model for calculating heat conduction losses in laser cutting processes. The key result is a correlation for the dimensionless rate of conduction heat loss into the workpiece: Q*cond = 3.20Pe+0.868, valid for the Peclet number range 0.20 ⱕ Pe ⱕ 10. The Peclet number is the primary dimensionless thermophysical parameter characterizing the cutting process, and may be considered to be a dimensionless cutting speed. Modern cutting speeds typically correspond to the range 0.3 ⱕ Pe ⱕ 1.5. In this range, the conduction loss increases with Pe more rapidly than does the enthalpy flux needed for melting the workpiece, which is exactly linear in Pe. At higher values of Pe, the opposite trend occurs. Additional correlations for characterizing the size of the HAZ show it to be extremely narrow in the immediate neighbourhood of the erosion front, typically on the order of 0.1–1.0 mm. The model was used to provide detailed temperature field results in this region, giving a clear picture of the interior of the HAZ. Farther from the erosion front, the model predictions for the temperature field are shown to be in good agreement with thermocouple measurements. The correlation for conduction loss was incorporated into an erosion front energy balance to predict the cutting speed. In dimensionless form, this balance results in an implicit function for the cutting speed, Pe, as a function of AintQ*las, Ste, and ⌫. These additional dimensionless parameters are the absorbed laser power, Stefan number, Ste (the ratio of sensible to latent heating required to melt the workpiece); and the dimensionless heat of combustion, ⌫ (the ratio of heat of combustion to sensible heat). The values of Ste and ⌫ are set only by the thermophysical properties of the workpiece. Thus for a given material, dimensionless cutting speed is a function only of absorbed laser power. By correlating kerf thickness and integrated absorptivity with workpiece thickness, the dimensionless result was recast into dimensional form, yielding uo = uo(Qlas/d). The predicted speeds closely matched experiments to within the experimental precision (about 10%). Development of a correlation for integrated absorptivity, Aint, as a function of the erosion front angle, ␥, required a judicious combination of analysis and experimental measurements. More generally, the local absorptivity depends upon the index of refraction, n*, and other parameters specifying the shape of the erosion front (assumed to be an oblique, semi-cylindrical prism in

J.M. Prusa et al. / International Journal of Machine Tools & Manufacture 39 (1999) 431–458

457

this study). The absorptivity analysis revealed several interesting characteristics. First, maximum absorptivity occurs when ␥ is close to the Brewster angle ( ⬇ 86.0° in this study). Second, the erosion front angle decreases as the workpiece thickness decreases. Finally, the location of the beam axis relative to the erosion front has a profound impact on Aint. A perturbation that causes a displacement of only 0.1 mm from the nominal beam/erosion front location can lower the integrated absorptivity by a factor of two. This is sufficient to destroy what otherwise would have been a high quality cut.

Acknowledgements The authors would like to thank the National Science Foundation for financial support of this work through the grant NSF-DMI-9407735; and J.L. Stanford for helpful comments on the absorptivity computations.

References [1] C.K. Lim, Numerical modeling of reactive gas assisted laser cutting of metals, Ph.D. Dissertation, Iowa State University, 1995, 126p. [2] J. Powell, Laser Cutting, Springer, Berlin, 1993, pp. 215–225. [3] M. Vicanek, C. Simon, Momentum and heat transfer of an inert gas jet to the melt in laser cutting, J. Phys. D: Appl Phys. 20 (1987) 1191–1196. [4] W.M. Steen, Laser Material Processing, Springer, Berlin, 1991, pp. 75–76. [5] A.R. Moss, J.A. Sheward, The arc plasma cutting of non-metallic materials, in: I.E.E.E. Conference on Electrical Methods of Machining, Forming, and Coating, London, 1970, pp. 41–47. [6] A.A. Wells, Conduction heat transfer, Welding J. 31 (Suppl.) (1952) 263–266. [7] W.M. Steen, Laser Material Processing, Springer, Berlin, 1991, pp. 155–161. [8] D. Rosenthal, Mathematical theory of heat distribution during welding and cutting, Welding J. 20 (1972) 2205–2345. [9] J.N. Gonsalves, W.W. Duley, Cutting thin metal sheets with CW laser, J. Appl. Phys. 43 (1972) 4684–4687. [10] K.A. Bunting, G. Cornfield, Toward a general theory of cutting: a relationship between the incident power density and the cut speed, J. Heat Transfer 97 (1975) 116–121. [11] M.J. Hsu, Analytical and experimental studies of advanced laser cutting techniques, Ph.D. Dissertation, Iowa State University, 1992, p. 137. [12] W. Schulz, D. Becker, J. Franke, R. Kemmerling, G. Herziger, Heat conduction losses in laser cutting of metals, J. Phys. D: Appl. Phys. 26 (1993) 1357–1363. [13] R.A. Danforth, Carbon dioxide laser cutting: beam mode and oxygen purity effects and validation of conduction heat losses. M.S. Thesis, Iowa State University, 1997. [14] G. Venkitachalam, P.A. Molian, An experimental study of laser cutting of a thick steel plate using a line beam, Trans. NAMRI/SME 24 (1996) 265–270. [15] T.R. Goodman, Application of integral methods to transient nonlinear heat transfer, Adv. Heat Transfer 1 (1964) 85–96. [16] V.S. Arpaci, Conduction Heat Transfer, Addison-Wesley, New York, 1966. [17] L.S. Yao, J. Prusa, Melting and freezing, Adv. Heat Transfer 19 (1989) 1–95. [18] D. Petring, P. Abels, E. Beyer, Absorption distribution on idealized cutting front geometries, Proc. SPIE 1020 (1988) 123–131. [19] J.B. Walsh, Electromagnetic Theory and Engineering Applications, Ronald Press, New York, 1960.

458

J.M. Prusa et al. / International Journal of Machine Tools & Manufacture 39 (1999) 431–458

[20] P.C. Clemmow, J.P. Dougherty, Electrodynamics of Particles, Addison-Wesley, New York, 1969. [21] Y.S. Touloukian, R.W. Powell, C.Y. Ho, P.G.K. Klemens, Thermal Conductivity Metallic Elements and Alloys, 1, Thermophysical Properties of Matter, The TPRC Data Series, IFI/Plenum, New York, 1981. [22] Y.S. Touloukian, E.H. Buyco, Specific Heat Metallic Elements and Alloys, 4, Thermophysical Properties of Matter, The TPRC Data Series, IFI/Plenum, New York, 1981. [23] P. Di Pietro, Y.L. Yao, A numerical investigation into cutting front mobility in laser cutting, Int. J Mach. Tools Manufact. 35 (1995) 673–688. [24] M.J. Hsu, P.A. Molian, Thermochemical modeling in CO2 laser cutting of carbon steel, J. Mat. Science 29 (1994) 5607–5611.