- Email: [email protected]
Finite element modelling of the laser cutting process
Compurers & Srrucfwes Vol. 49, No. 2. pp. 231-241. 1993 Prinlcd in Great Britain.
FINITE
0
lws-7949193 16.00 + 0.00 1993 Rrgamon Press Ltd
ELEMENT MODELLING OF THE LASER CUTTING PROCESS
M. J. KIM, Z. H. CI-IENand P. MAJIJMDAR Department of Mechanical Engineering, Northern Illinois University, Dekalb, IL 60115, U.S.A. (Received 10 August 1992) Abstract-A
two-dimensional transient model for the removal of metal using a high power laser is developed. The nondimensional forms of the governing equation and the boundary conditions are derived, and the numerical solution is obtained by using a finite element method. This finite element model is used to study the performance of material removing process for a wide range of parameters. The performance characteristics and the effects of these parameters on the shape and size of the groove formed are discussed.
INTRODUCTION
studied by Soodak [l]. He considered constant heat-
of the most promising applications of high-power continuous-wave gas lasers is the cutting of materials. The cutting process is carried out by melting and/or evaporation of a material with a focused Gaussian laser beam as a heat source. The metal can be. cut by the total removal of the material along the cutting line or by partial removal by processes such as scribing, or cleaving used particularly for brittle materials. A sketch of a typical laser cutting installation is shown in Fig. 1. The laser beam is focused by a lens through which the laser beam passes and behind which it converges to a minimum beam waist around the focal point of the lens. A jet of gas is directed parallel to the laser beam on the material. The gas is supplied through a conical nozzle. The laser beam also passes through the nozzle. With this concentrated energy at the focal point of the lens, it is possible to heat, melt and vaporize any material. The laser beam produces a heat source with a high energy concentration; therefore, it provides a smaller cutting width, smaller sizes of the heat affected zone, smaller mechanical effects in the material, and higher cutting rate than any other thermal cutting techniques. The width of the cutting zone is close to the diameter of the laser-focused spot. Thus this technique offers a high tolerance, rapid, flexible and precise cutting tool which is non-contacting and is not affected by the hardness of the materials. However, it is very important to the success of this process that the performance, quality, and thermal effects of the process is well understood. Most of the theoretical works on laser-treatment heat transfer to data have centered on the solution of the classical heat conduction equation for a stationary or moving semi-infinite solid. Cases with and without phase change, and a variety of irradiation or source conditions have been studied. The problem of melting with complete removal of melt subjected to different types of boundary conditions has been
ing of the surface and numerically evaluated the steady state melting rate. Landau [2] considered melting with complete removal of the melt from a one-dimensional slab with one end subjected to time varying heating and the other end insulated. He obtained the time-dependent temperature distribution in the solid, the position of the melting surface, and the limiting speed of the melting front. Rogerson and Chayt [3] calculated the total melting time of a one-dimensional ablating slab subjected to constant heating with complete removal of melt. They showed that the processing is independent of the transport properties of the material. Dabby and Paek [4] have analytically studied the problem of material removal from the surface by considering laser penetration into the solid that is vaporizing on its free surface. They also assumed that the gas created by the vaporization is transparent to the incident laser beam. Von Allmen [5] considered drilling processes with material expulsion due to pressure gradients caused by evaporation. His theoretical results showed considerable agreement with experiments. Cline and Anthony [6] derived a model for laser heating and melting of materials for a Gaussian source moving at constant velocity, for which they calculated the temperature distribution and depth of the melting zone as a function of laser beam diameter, velocity and power. A three-dimensional heat transfer model was developed by Mazumder and Steen [7] for a laser beam striking the surface of an opaque substrate moving with uniform velocity. The model was solved by finite difference method and results were presented for temperature distribution and melt depth. Modest and Abakians [8] studied numerically the formation of a groove by evaporation on a moving semi-infinite solid. Results for groove depth and shape were presented for a variety of laser and solid parameters. Hsu et al. [9] investigated the effect of a continuous wave laser on surface layer melting and
One
231
232
M. J. KIM er al. ~ianla8erbeam
Pocusillg lens
the intensities required for laser drilling, and the purpose of this model was to determine the effect of spatially varying laser intensities on the temperature distribution. The model accounts for both sensible heating and phase changes (solid-to-liquid and liquidto-vapor). Glass et al. [14] investigated the effects of various parameters in laser cutting of metallic glass ribbons which undergo ductile-to-brittle phase transitions when heated above crystallization temperatures. They modeled the laser/material interaction using a quasi-steady, three-dimensional finite difference technique to predict the temperatures and cooling rates in the heat-affected zone and compared this with experimental results. The objective of this work is to develop a finite element model for a material removal process using a high energy laser, and to use this model to carry out a parametric study.
Fig. 1. Laser cutting installation. MATHEMATICAL
subsequent solidification using a one-dimensional heat conduction model. His model calculates interface velocities of melting and solidification in addition to the calculation of temperature profiles and rates of heat transfer. Biyikli and Modest [lo] studied the effect of focusing parameters on the formation of groove depth and shape using a numerical model. Chan et al. [ 1l] developed one-dimensional transient and steady state models describing the process of material removal by vaporization and liquid expulsion. The problem of moving the solid-liquid interface and the vapor-liquid interface are modeled by boundary immobilitation and a Mot&-Smith type solution. Grigoropoulos et al. [ 121 described detailed experimental observations of the associated phase change process and presented a computational conductive heat transfer model for laser melting and recrystallization of thin semiconductor films. Minardi and Bishop [ 131 developed a two-dimensional transient computer model to determine the temperature distribution within a material subjected to irradiation for
FORMULATION
The mathematical model describing the process of material removal from the surface subjected to high intensity laser beam is presented here. Figure 2 shows the schematic diagram of the process of material removal with different regimes. Laser and optical parameters
For any expanding laser beam striking a surface at any arbitrary angle, the beam intensity is expressed as F(x, y, z) = (L + tan 6 cos f#~i^ + tan 0 sin +i)F,
Ri
-(x2 + ,WRZ(r)
xR2(2)Xe
where i,j, k are unit vectors in x, y and z directions, 0 is the angle between the direction of laser beam and z axis, #I is the angle made by the laser beam with the x axis, F, = (P/xRi) is the beam power density at beam center, P is the laser power, & is the effective laser beam radius at the focal plane, and R(0) is the
Fig. 2. Schematic diagram of three regimes.
Finite element modelling of the laser cutting
beam radius at the material surface, z = 0. The radius of the expanding laser beam is given as R(z)=&[l
+(%>,I,,,.
(2)
AtZ=S(X,
-cymax.
process
Y),x,,(Y)
a0 do = O or
0 =0
At Z=Z,,
Heat transfer model
Boundary condition 3
in the domain o
-co
+ Bi8 = 0
Y
(7)
Y< Y,.
0=0
(8)
At X = +X,. Boundary condition 4
ao
=+B,0=0
or
At Y= Y,, -X,
@=O
(9)
O,
Boundary condition 5
a8
az-BiO=O
AtZ=O,
(IO)
X,,,,(Y)
-Y,,,,
Boundary condition 6
az
-cc
g
-Yr<
as-B.0=0
(3)
or
-X,
Governing equations
p+“~=~+E?+~
-Y,,,,<
Boundary condition 2
where I is the wave length of laser beam and W is the distance between focal plane of concentrating optical lens and the material surface.
The process of material removal by high intensity laser can be described by a heat transfer model based on the conservation of energy. It is assumed that there are two different regions on the surface of the material: One is assumed to be the area in which evaporation takes place and the other is too far away to have reached evaporation temperature. The following assumptions are also made in deriving the heat transfer model: (1) material moves at constant relative velocity, (2) solid is isentropic with constant thermal properties, (3) material is opaque and has constant absorptivity, (4) evaporated material does not interfer with the incident laser beam, (5) phase change from solid to vapor occurs in one step, (6) convection and radiation losses from the surfaces to environment can be approximated by using a single constant coefficient. By using these assumptions, the conservation of energy over an infinitesimal element gives the mathematical statement of the problem as follows:
233
At Y= -Yr,
or
1
@=O
(11)
-X,
Dimensionless variables
BOUNDARY CONDITIONS
Boundary condition 1
Z
UT
(12)
Z=R(0)’f=RZ(o).
Region I
Parameters iVk(BiO -g)
= e--(x2+y2).
(4) u
NO)
=
-9
At Z=O,
-X,
-yF<
YZF yr=R(0)
Y< Y,. zF=;,
Region II
X=XZF F R(0)’
N
k
B+!?
=W,,,-
Tm)
R(O)aF,
’
phcu Ne= aF, ’
(13)
FINITE ELEMENT FORMULATION
In order to develop a finite element formulation 0 =I.
(6) of the problem, we assumed the following general
234
M.
J. KIM et al.
differential equation and boundary conditions governing the transient heat transfer in a twodimensional region R and total boundary r as
a8
Equation
(20) can be expressed in matrix form as [M’,“]{0} + [X$‘]{0} = {FI”},
(21)
where
dt+U~-&(;)--$)+/=O
(14)
inR,O
a0 Enx+~yny-BiO +q=O
on r,
(1%
and 0 = 0,
on r2.
(16)
Initial conditions 0=0,innwhent=O
(17)
where Bi, @,, f and q are given functions of position and time. The variational formulation of eqn (14) over an element ace) is obtained by multiplying eqn (14) by test function v and integrating.
a0
vat+Uv-+
-
a0
a8 av
a8ay
-
I
BiYiYj& r,r,
fi’=
-
Q(*, Y)=C,+C*X+C,Y,
where C, , C,, and C, are known functions of nodal values 0i, 0r, 0, and global coordinates of the three vertices of the triangle. Rearranging eqn (22) we may express 0 as a function of nodal values and global coordinates as
a0
QGK Y) = Q, Yu,(-KY) + Q,Yu,(X Y)
a8
ax’““+arvn’ > SC
ds =o. (18)
0
can be approximated
C
j-1
by the
yh
0j(r)Yj(X9
-
g,Y,z s MWI)
8, u (
dXdY
where the linear interpolation triangular element is given by
(19)
Y I”) =
av,ayi
j-l
O(r)
J axax+arar
O,BiYiYj& s IW,
&
8iyyI”,
(23)
function
YUp)for a
(ai + pi* + yi Y)
(24)
e Xj Y,- X, Y,
(25)
f3i =
r, - Y,
(26)
yi =
x, - xj
(27)
ai =
2A,=
1
Y)= i
+0jY,(x,
where Sj are the nodal values of 0 and Y,(X, Y) are the interpolation functions. Substituting eqn (19) for 0 and v = Yi into eqn (18), we get
+i
(22)
ax
s(x~ t, =
+u
s r(r)
When three node triangular elements are used, the unknown function 0 may be approximated as varying linearly within the element as
r(c)
The function following form
qYi ds.
Y,fdXdY-
s IWC)
1 Xi Y, 1 X,Y, 1 X, Y, (
=(X*X,-X,Y,)+(X,Y,-X,Y,)
ayiauj d*dY
+(X,Y,-XZY,).
1
In order to compute the coefficients of eqn (21), those coefficients are rewritten as [MI;‘1 = ]%I
+~~~~~qY,~+I_,Y,fdXdY=O.
(28)
(20) @‘I = %]&I + a,,Nji]
+ a,2[S$? + a&;*]
(29) (30)
235
Finite element modelling of the laser cutting process
{fl’) = [#I + IQ?%
(31)
where ‘P,!P,dX dY,
S,=
s:/’ =
s
!!?alv,dXdY
a*,
S’2 r/ =
ax dX
LI!Piz s oW
[A]{&}+[B]{@}={P}
(32)
dX dY,
(33)
fl) =
Using interelement continuity conditions and correspondence between the local nodes and the global nodes given by boolean connectivity matrix into eqn (41), we get the global finite element model of the problem in matrix as
j-J?Pi dX d Y
(34)
qnYi ds. s I-W
(35)
Using the linear interpolation function given by eqn (24) in eqns (32)-(34) we obtain S+
Sb’
(42)
where [A I, PI, {P} are known matrices and (0) is the column vector of the unknowns. The set of first order differential equations can now be transformed into a set of algebraic equations by approximating the first order time derivative. The first order time derivatives can be approximated by using a scheme given as
J@L+rw” Qy’ =
forO
At,+,
for 0 < c < 1.
(43)
Using this approximation scheme, eqn (42) is transformed into a set of algebraic equations given in matrix form as
(36) where
[B]= [A]- (I- c)At,+,[Bl. (37) + f Cr?JJ8iSj+l*l(YtSj + YjSt)+ 1Ci2YiY/ll
f@~=J(a,+fl~+y.P)=fA (38) I I I
2
3’
For the boundary integral given by the eqn (35), it is only necessary to evaluate the integral over the portion of r’ that falls on the boundary using a one dimensional element with a one-dimensional interpolation function as h qnyih
s
Q,=
q,oh
9
(39)
i 2 where the interpolation Y ,=1-s/h,
function is given by Iyz=s/n.
(40)
The assembly of all the elements is obtained by assuming that the quadratic functional associated with the problem is the summation of the quadratic functionals of elements as below
The solution at time t = t, + , is obtained in terms of the solution known at time t, by inverting the matrix [A]. At t = 0 the solution is known from the initial conditions of the problem, and therefore, eqn (44) can be used to obtain the solutions with marching in time. COMPUTATIONALMETHODS
The problem domain is initially divided up into a mesh of triangular finite elements as illustrated in Fig. 3. Since the center regions will be subjected to higher temperature gradients as they are located right below the laser beam, smaller elements are used in this section. However, as the location of the surface at 2 = S moves due to the evaporation of the materials, we need to relocate this surface based on the criteria that evaporation takes place when temperature at any node reaches unity, and redistribute the meshes in the new domain. In order to establish the domain at every time step an iterative procedure, as illustrated in Fig. 4, is utilized. Figure 4 shows the iterative scheme to determine the interface location S. At any section a new position of interface location Z,,, = S is identified and calculated based on linear interpolation of two
M. J. KIMer al.
236
I
3
2
4
5 6 7 6 9 1011121314l51617
16
19 R
20
21 1
Fig. 3. Two-dimensional finite element mesh. other points between which the value of 8 crosses over. This position is expressed as below zpoin,= z3 + z,-s *4 _ & (8, - 1).
--X
--
____--___---
+-IY4
y*
Y4
1
Y
e,=i.i
y3
p
1.3
e,=
y2
f
1.4
e,=
Yl
____
(45)
Once Zpoinris identified, the length between 2 = S and ZE) is divided into equal lengths among the elements and then this new mesh distribution is used to obtain a new temperature distribution. This iteration procedure will be continued until a convergence
Y*
I
---p-
e, = 0.9
Ql
--
y3
-
_
* \
-
-I-
I
‘1
I
‘1
X
e4
93
6
j--+--t Y
Fig. 4. Iteration scheme.
I
1. t = 0.3 2. t = 0.4 3. t = 0.5 4. t = 0.6 5. 6. tt = = 0.7 0.8
-o.sl 4.0 -2.5 I
”
I
I
I
I
I
I
I
I
I
-1.0
-0.7
-0.3
0
0.3
0.7
1.0
2.5
4.0
Porition along the
length x
Fig. 5. The depth of the groove at different time t (LI = 0, N, = 0, Bi = 0.0001, Nk = 0.5).
PTP
237
Finite element modelling of the laser cutting process
-1.4 I -4.0
I -2.5
I -1.0
I I I I I -0.1 -0.3 0 0.3 0.7 Position along the length X
I 1.0
I 2.5
I 4.0
Fig. 6. The depth of the groove at different time I (CI = 0, N, = 0, Bi = 0.0001, Nk = 0.5).
in temperature at the surface is achieved within a certain assigned accuracy limit. RESULTS AND DISCUSSIONS
The finite element model is used to study the process of material removal from a solid surface subjected Gaussian laser beam of constant power. The evaporation rate, groove shape and temperature distribution in the media are studied for both steady state and transient conditions. Solutions are obtained for various operating parameters and cases with different boundary conditions, Transient state results Transient phenomena of material removal were studied first. Figures 5 and 6 show the variation of
depth of the groove along the axial position x with an increase in time for the case of U = 0. However, the rate of increment decreases as time increases and the heat losses are increased with an increase in area of interface. Figure 6 shows the final depth of the groove which remains steady with an increase in time. Figure 7 shows the variation of groove depth with an increase in time for the case with U = 1. The results show that the groove shape is similar to that in the case of U = 0, except it is slightly asymmetrical and reaches a stable position faster due to smaller residence time. Figures 8 and 9 show the comparison of groove shape given by both steady state and transient state analysis for the cases U = 0 and U = 1, respectively. It can be seen that this unsteady state solution is very close to the solution given by steady
m -0.6 t
-0.7-0.8 -“*9 -1.0 -1.1 -4.0
1.tm0.4 2. t = 0.5 3. t = 0.6 4. t = 0.7 5. t = 0.8 6. t = 4.5 I -2.5
I
I
-1.0
-0.7
6
,Vl -0.3
0
0.3
I 0.7
1.0
2.5
4.0
Position along the length x Fig. 7. The depth of the groove at different time r (U = I, N. = 0.01, Ei = 0.0001, Nk = 0.5). CAS 49,2-c
0 -0.1 -0.2 -0.3 -0.4 -0.5 -0.6 3
-0.7
*
-0.8 -0.9 -1.0 -1.1 -1.2
1.
unsteady state: t = 7.9
-1.3 -1.4 -4.0
-2.5
-1.0
-0.7
-0.3
0
0.3
0.7
1.0
2.5
10
Position along the length X
Fig. 8. The depth of the groove for steady state and unsteady state for U = 0, N, = 0, Bi = 0.0001, Nk = 0.5.
-0.1 -0.2 -0.3 -0.4 R
-0.5 -
:
-0.6 -0.7 -0.8 -0.9 -1.0 -1.1 -4.0
1. unsteady state: t =
4.5
2. steady state
I -2.5
I -1.0
I
I
I
I
I
I
-0.7
-0.3
0
0.3
0.7
1.0
Position along the
2.5
4.0
length X
Fig. 9. The depth of the groove for steady state and unsteady state for CJ= 1, N, = 0.01, Bi = 0.0001. Nk = 0.5.
B.C.:
B=Oatx=fq
Position along the length x Fig. 10. The depth of the groove for the different values of Nk. 238
Finite clement modelling of the laser cutting process
239
-0.8 -0.9 -
-1.0 _t.t _ -1.2 -
1. Bi 2. Bi 3. Bi 4. Bi
5 0.0001 re:0.01 = 0.1 = 0.3
5. Bi - 0.7 6. Bi = 0.9
-1.3 -4.0
-2.5
-1.0
-0.1
-0.3
0
0.3
0.7
1.0
2.5
4.0
Position along the length X Fig. 11. The depth of the groove for the different values of
state ana!ysis. The above results validate the usefulness of the steady state solutions. Therefore, the material removal rate and given shape can be predicted with good approximation by using steady state analysis.
Steady state results The steady state model is further used to study the material removal process with varying parameters and boundary conditions. Figures IO-16 show the results assuming the boundary condition at 2 = Z,, as 8 = 0. Figure 10 shows the variation of groove shapes with the different values of Nk. It can be seen from the definition of N,, that it is the ratio of conduction heat loss along the interface to the absorbed laser power. For large values of N,, the
Bi.
laser energy is mostly taken away by conductive heat losses and this causes the formation of the shallow groove. Decreasing Nk will cause a reduction in conductive heat losses and increases in the depth of the groove. However, as the groove becomes deeper, the radiation per unit surface area becomes smaller, resulting in less energy transfer into the material. At this stage, increasing Nk will not give higher material removal rate. In order to see the effects of surface convective heat losses on the groove shape and depth, results are presented in Fig. 11 for Nk = 0.5 with varying Biot numbers. For a fixed conduction resistance, a decrease in Biot number corresponds to increase in surface convection resistance, and this causes higher material removal rate. It can be seen from the results
B.C.=g=Oatr=fao,
-4.0
-2.3
-1.0
-0.7
-0.3
0
0.3
0.7
1.0
2.5
Porftion along the length x
Fig. 12. The depth of the groove for the different values of U.
1
-0.2 -0.4 ,+e+-+-+-+
-0.6 -0.8 -1.0 h -1.2 5 N -1.4 -1.6 -1.8 -2.0
x-x-x-.x-~
_w~xcx-x-x-x-x
-2.2
-1.0 -0.7
-0.3
0
0.3
0.7
1.0
2.5
4.0
Position along the length X Fig. 13. The distribution of the node points along the x-direction for U - 0, N, = 0, Nk = 0.5, Bi = 0.0001.
Position along the length X Fig. 14. The temperature distribution along the X direction for .!I = 0, N, = 0, Nk = 0.5, Bi = 0.0001. 0
-0.2 -0.4 -0.6 -0.8 -1.0 fi -1.2 5 N -1.4 -1.6 -1.8 -2.0
x-x-x-%-x_y~
_-x-x-x-x
-2.2
-0.7
-0.3
0
0.3
0.7
1.0
2.5
4.0
Position along the length x Fig. IS. The distribution of the node points along the X direction for U = 1, N, = 0.01, Nk = 0.5, Bi = 0.0001. 240
241
Finite element modelling of the laser cutting process
-1.0
-0.7
-0.3
0
0.3
0.7
1.0
1.5
4.0
Position along the length X 16. The temperature distribution along the X direction for U = 1, N, = 0.01, Nk = 0.5, Bi = 0.0001.
that groove depth and shape increase with a decrease in Biot number. Results also show that no appreciable change occurs in the groove if the Biot number is less than 0.01. Figure 12 shows the effect of nondimensional velocity U on the shape and material removal rate, shape of groove will change with the change of nondimensional velocity U, it relates the relative speed of the laser to that of thermal diffusion into the medium. For larger values of U material removal rate will be less because of small residence time. Therefore, both the groove depth and the slope of the left side edge will be decreased with an increase in value of U. Figure 13 shows the shape of groove and the positions of the finite element node points for the case with U = 0. As the material evaporates during the cutting process, meshes are regenerated in the media and nodes along each layer are displaced from their original position. Figure 14 shows the nodal temperature distribution along each layer. It can be seen that at the top layer the temperature of all nodes approach the evaporation temperature and this is uniform along throughout. At other layers, the nodes located away from the center of laser focal point have a lower temperature than those located at the center of the focal region. This is because less heat is taken away by conduction in the media than that used in evaporation for the material at the center node. Figures 15 and 16 show similar results for the case with U = 1. The groove shape and the temperature distribution is no longer symmetrical, and also less material is removed due to lower residence time. REFERENCES I. H. Soodak, Effects of heat transfer between gases and solids. Ph.D. thesis, Duke University, Durham, NC (1943).
2. H. G. Landau, Heat conduction in a melting solid. Q. appl Math. 8, 81-94 (1950). 3. J. E. Rogerson and G. A. Chayt, Total melting time in the ablating slab problem. J. appl. Phys. 42, 27 1I-27 13 (1971). 4. F. W. Dabby and U.-C. Paek, High-intensity laser-induced vaporization and explosion of solid material. IEEE J. Quantum Electronics QE-8, 106-l 11 (1972).
5. M. Von Allmen, Laser drilling velocity in metajs. J. appl. Phys. 47, 5460-5463
(1976).
6. H. E. Cline and T. R. Anthony, Heat treating and melting material with scanning laser or electron beam. J. appl. Phys. 48, 3895-3900 (1977). 7. J. Mazumder and W. M. Steen, Heat transfer model for CW laser material processing. J. appl. Phys. 51,941-947 (1980). 8. M. F. Modest and H. Abakians, Evaporative cutting of a semi-infinite body with a moving CW laser. Trans. ASME 108. 602-607 (1986). 9. S. C. Hsu, S. Chakravorty aid R: Mehrabian, Rapid melting and solidification of a surface layer. Metall. Trans.
9B, June (1978).
10. S. Biyikli and M. F. Modest, Effect of beam expansion
and focusing on evaporative cutting with a moving CW laser. J. Heal Transfer 110, 529-532 (1988). 11. C. Chan, J. Mazumdar and M. M. Chen, A two dimensional transient model for convection in laser melted pool. Mefall. Trans. lSA, 2175-2184 (1984). 12. C. P. Grigoropoulos, S. E. Long, A. F. Emery and W. E. D&her, Jr, Experimental and computational analysis of laser melting of thin silicon films. 1989 National Heat Transfer Conference, HTD-Vol. 113, Heat Transfer in Manufacturing and Materials Processing. 13. A. Minardi and P. J. Bishop, Temperature distribution within a metal subjected to irradiation by a laser of spatially varying intensity. 1989 National Heat Transfer Conference, HTD-Vol. 113, Heat Transfer in Manufacturing and Materials Processing. 14. J. M. Glass, H. P. Groger, R. J. Churchill, J. W. Lindau and T. E. Diller, Heat transfer in metallic glasses during laser cutting. 1989 National Heat Transfe; Conference, HTD-Vol. 113, Heat Transfer in Manufacturing and Materials Processing.
FINITE
0
lws-7949193 16.00 + 0.00 1993 Rrgamon Press Ltd
ELEMENT MODELLING OF THE LASER CUTTING PROCESS
M. J. KIM, Z. H. CI-IENand P. MAJIJMDAR Department of Mechanical Engineering, Northern Illinois University, Dekalb, IL 60115, U.S.A. (Received 10 August 1992) Abstract-A
two-dimensional transient model for the removal of metal using a high power laser is developed. The nondimensional forms of the governing equation and the boundary conditions are derived, and the numerical solution is obtained by using a finite element method. This finite element model is used to study the performance of material removing process for a wide range of parameters. The performance characteristics and the effects of these parameters on the shape and size of the groove formed are discussed.
INTRODUCTION
studied by Soodak [l]. He considered constant heat-
of the most promising applications of high-power continuous-wave gas lasers is the cutting of materials. The cutting process is carried out by melting and/or evaporation of a material with a focused Gaussian laser beam as a heat source. The metal can be. cut by the total removal of the material along the cutting line or by partial removal by processes such as scribing, or cleaving used particularly for brittle materials. A sketch of a typical laser cutting installation is shown in Fig. 1. The laser beam is focused by a lens through which the laser beam passes and behind which it converges to a minimum beam waist around the focal point of the lens. A jet of gas is directed parallel to the laser beam on the material. The gas is supplied through a conical nozzle. The laser beam also passes through the nozzle. With this concentrated energy at the focal point of the lens, it is possible to heat, melt and vaporize any material. The laser beam produces a heat source with a high energy concentration; therefore, it provides a smaller cutting width, smaller sizes of the heat affected zone, smaller mechanical effects in the material, and higher cutting rate than any other thermal cutting techniques. The width of the cutting zone is close to the diameter of the laser-focused spot. Thus this technique offers a high tolerance, rapid, flexible and precise cutting tool which is non-contacting and is not affected by the hardness of the materials. However, it is very important to the success of this process that the performance, quality, and thermal effects of the process is well understood. Most of the theoretical works on laser-treatment heat transfer to data have centered on the solution of the classical heat conduction equation for a stationary or moving semi-infinite solid. Cases with and without phase change, and a variety of irradiation or source conditions have been studied. The problem of melting with complete removal of melt subjected to different types of boundary conditions has been
ing of the surface and numerically evaluated the steady state melting rate. Landau [2] considered melting with complete removal of the melt from a one-dimensional slab with one end subjected to time varying heating and the other end insulated. He obtained the time-dependent temperature distribution in the solid, the position of the melting surface, and the limiting speed of the melting front. Rogerson and Chayt [3] calculated the total melting time of a one-dimensional ablating slab subjected to constant heating with complete removal of melt. They showed that the processing is independent of the transport properties of the material. Dabby and Paek [4] have analytically studied the problem of material removal from the surface by considering laser penetration into the solid that is vaporizing on its free surface. They also assumed that the gas created by the vaporization is transparent to the incident laser beam. Von Allmen [5] considered drilling processes with material expulsion due to pressure gradients caused by evaporation. His theoretical results showed considerable agreement with experiments. Cline and Anthony [6] derived a model for laser heating and melting of materials for a Gaussian source moving at constant velocity, for which they calculated the temperature distribution and depth of the melting zone as a function of laser beam diameter, velocity and power. A three-dimensional heat transfer model was developed by Mazumder and Steen [7] for a laser beam striking the surface of an opaque substrate moving with uniform velocity. The model was solved by finite difference method and results were presented for temperature distribution and melt depth. Modest and Abakians [8] studied numerically the formation of a groove by evaporation on a moving semi-infinite solid. Results for groove depth and shape were presented for a variety of laser and solid parameters. Hsu et al. [9] investigated the effect of a continuous wave laser on surface layer melting and
One
231
232
M. J. KIM er al. ~ianla8erbeam
Pocusillg lens
the intensities required for laser drilling, and the purpose of this model was to determine the effect of spatially varying laser intensities on the temperature distribution. The model accounts for both sensible heating and phase changes (solid-to-liquid and liquidto-vapor). Glass et al. [14] investigated the effects of various parameters in laser cutting of metallic glass ribbons which undergo ductile-to-brittle phase transitions when heated above crystallization temperatures. They modeled the laser/material interaction using a quasi-steady, three-dimensional finite difference technique to predict the temperatures and cooling rates in the heat-affected zone and compared this with experimental results. The objective of this work is to develop a finite element model for a material removal process using a high energy laser, and to use this model to carry out a parametric study.
Fig. 1. Laser cutting installation. MATHEMATICAL
subsequent solidification using a one-dimensional heat conduction model. His model calculates interface velocities of melting and solidification in addition to the calculation of temperature profiles and rates of heat transfer. Biyikli and Modest [lo] studied the effect of focusing parameters on the formation of groove depth and shape using a numerical model. Chan et al. [ 1l] developed one-dimensional transient and steady state models describing the process of material removal by vaporization and liquid expulsion. The problem of moving the solid-liquid interface and the vapor-liquid interface are modeled by boundary immobilitation and a Mot&-Smith type solution. Grigoropoulos et al. [ 121 described detailed experimental observations of the associated phase change process and presented a computational conductive heat transfer model for laser melting and recrystallization of thin semiconductor films. Minardi and Bishop [ 131 developed a two-dimensional transient computer model to determine the temperature distribution within a material subjected to irradiation for
FORMULATION
The mathematical model describing the process of material removal from the surface subjected to high intensity laser beam is presented here. Figure 2 shows the schematic diagram of the process of material removal with different regimes. Laser and optical parameters
For any expanding laser beam striking a surface at any arbitrary angle, the beam intensity is expressed as F(x, y, z) = (L + tan 6 cos f#~i^ + tan 0 sin +i)F,
Ri
-(x2 + ,WRZ(r)
xR2(2)Xe
where i,j, k are unit vectors in x, y and z directions, 0 is the angle between the direction of laser beam and z axis, #I is the angle made by the laser beam with the x axis, F, = (P/xRi) is the beam power density at beam center, P is the laser power, & is the effective laser beam radius at the focal plane, and R(0) is the
Fig. 2. Schematic diagram of three regimes.
Finite element modelling of the laser cutting
beam radius at the material surface, z = 0. The radius of the expanding laser beam is given as R(z)=&[l
+(%>,I,,,.
(2)
AtZ=S(X,
-cymax.
process
Y),x,,(Y)
a0 do = O or
0 =0
At Z=Z,,
Heat transfer model
Boundary condition 3
in the domain o
-co
+ Bi8 = 0
Y
(7)
Y< Y,.
0=0
(8)
At X = +X,. Boundary condition 4
ao
=+B,0=0
or
At Y= Y,, -X,
@=O
(9)
O,
Boundary condition 5
a8
az-BiO=O
AtZ=O,
(IO)
X,,,,(Y)
-Y,,,,
Boundary condition 6
az
-cc
g
-Yr<
as-B.0=0
(3)
or
-X,
Governing equations
p+“~=~+E?+~
-Y,,,,<
Boundary condition 2
where I is the wave length of laser beam and W is the distance between focal plane of concentrating optical lens and the material surface.
The process of material removal by high intensity laser can be described by a heat transfer model based on the conservation of energy. It is assumed that there are two different regions on the surface of the material: One is assumed to be the area in which evaporation takes place and the other is too far away to have reached evaporation temperature. The following assumptions are also made in deriving the heat transfer model: (1) material moves at constant relative velocity, (2) solid is isentropic with constant thermal properties, (3) material is opaque and has constant absorptivity, (4) evaporated material does not interfer with the incident laser beam, (5) phase change from solid to vapor occurs in one step, (6) convection and radiation losses from the surfaces to environment can be approximated by using a single constant coefficient. By using these assumptions, the conservation of energy over an infinitesimal element gives the mathematical statement of the problem as follows:
233
At Y= -Yr,
or
1
@=O
(11)
-X,
Dimensionless variables
BOUNDARY CONDITIONS
Boundary condition 1
Z
UT
(12)
Z=R(0)’f=RZ(o).
Region I
Parameters iVk(BiO -g)
= e--(x2+y2).
(4) u
NO)
=
-9
At Z=O,
-X,
-yF<
YZF yr=R(0)
Y< Y,. zF=;,
Region II
X=XZF F R(0)’
N
k
B+!?
=W,,,-
Tm)
R(O)aF,
’
phcu Ne= aF, ’
(13)
FINITE ELEMENT FORMULATION
In order to develop a finite element formulation 0 =I.
(6) of the problem, we assumed the following general
234
M.
J. KIM et al.
differential equation and boundary conditions governing the transient heat transfer in a twodimensional region R and total boundary r as
a8
Equation
(20) can be expressed in matrix form as [M’,“]{0} + [X$‘]{0} = {FI”},
(21)
where
dt+U~-&(;)--$)+/=O
(14)
inR,O
a0 Enx+~yny-BiO +q=O
on r,
(1%
and 0 = 0,
on r2.
(16)
Initial conditions 0=0,innwhent=O
(17)
where Bi, @,, f and q are given functions of position and time. The variational formulation of eqn (14) over an element ace) is obtained by multiplying eqn (14) by test function v and integrating.
a0
vat+Uv-+
-
a0
a8 av
a8ay
-
I
BiYiYj& r,r,
fi’=
-
Q(*, Y)=C,+C*X+C,Y,
where C, , C,, and C, are known functions of nodal values 0i, 0r, 0, and global coordinates of the three vertices of the triangle. Rearranging eqn (22) we may express 0 as a function of nodal values and global coordinates as
a0
QGK Y) = Q, Yu,(-KY) + Q,Yu,(X Y)
a8
ax’““+arvn’ > SC
ds =o. (18)
0
can be approximated
C
j-1
by the
yh
0j(r)Yj(X9
-
g,Y,z s MWI)
8, u (
dXdY
where the linear interpolation triangular element is given by
(19)
Y I”) =
av,ayi
j-l
O(r)
J axax+arar
O,BiYiYj& s IW,
&
8iyyI”,
(23)
function
YUp)for a
(ai + pi* + yi Y)
(24)
e Xj Y,- X, Y,
(25)
f3i =
r, - Y,
(26)
yi =
x, - xj
(27)
ai =
2A,=
1
Y)= i
+0jY,(x,
where Sj are the nodal values of 0 and Y,(X, Y) are the interpolation functions. Substituting eqn (19) for 0 and v = Yi into eqn (18), we get
+i
(22)
ax
s(x~ t, =
+u
s r(r)
When three node triangular elements are used, the unknown function 0 may be approximated as varying linearly within the element as
r(c)
The function following form
qYi ds.
Y,fdXdY-
s IWC)
1 Xi Y, 1 X,Y, 1 X, Y, (
=(X*X,-X,Y,)+(X,Y,-X,Y,)
ayiauj d*dY
+(X,Y,-XZY,).
1
In order to compute the coefficients of eqn (21), those coefficients are rewritten as [MI;‘1 = ]%I
+~~~~~qY,~+I_,Y,fdXdY=O.
(28)
(20) @‘I = %]&I + a,,Nji]
+ a,2[S$? + a&;*]
(29) (30)
235
Finite element modelling of the laser cutting process
{fl’) = [#I + IQ?%
(31)
where ‘P,!P,dX dY,
S,=
s:/’ =
s
!!?alv,dXdY
a*,
S’2 r/ =
ax dX
LI!Piz s oW
[A]{&}+[B]{@}={P}
(32)
dX dY,
(33)
fl) =
Using interelement continuity conditions and correspondence between the local nodes and the global nodes given by boolean connectivity matrix into eqn (41), we get the global finite element model of the problem in matrix as
j-J?Pi dX d Y
(34)
qnYi ds. s I-W
(35)
Using the linear interpolation function given by eqn (24) in eqns (32)-(34) we obtain S+
Sb’
(42)
where [A I, PI, {P} are known matrices and (0) is the column vector of the unknowns. The set of first order differential equations can now be transformed into a set of algebraic equations by approximating the first order time derivative. The first order time derivatives can be approximated by using a scheme given as
J@L+rw” Qy’ =
forO
At,+,
for 0 < c < 1.
(43)
Using this approximation scheme, eqn (42) is transformed into a set of algebraic equations given in matrix form as
(36) where
[B]= [A]- (I- c)At,+,[Bl. (37) + f Cr?JJ8iSj+l*l(YtSj + YjSt)+ 1Ci2YiY/ll
f@~=J(a,+fl~+y.P)=fA (38) I I I
2
3’
For the boundary integral given by the eqn (35), it is only necessary to evaluate the integral over the portion of r’ that falls on the boundary using a one dimensional element with a one-dimensional interpolation function as h qnyih
s
Q,=
q,oh
9
(39)
i 2 where the interpolation Y ,=1-s/h,
function is given by Iyz=s/n.
(40)
The assembly of all the elements is obtained by assuming that the quadratic functional associated with the problem is the summation of the quadratic functionals of elements as below
The solution at time t = t, + , is obtained in terms of the solution known at time t, by inverting the matrix [A]. At t = 0 the solution is known from the initial conditions of the problem, and therefore, eqn (44) can be used to obtain the solutions with marching in time. COMPUTATIONALMETHODS
The problem domain is initially divided up into a mesh of triangular finite elements as illustrated in Fig. 3. Since the center regions will be subjected to higher temperature gradients as they are located right below the laser beam, smaller elements are used in this section. However, as the location of the surface at 2 = S moves due to the evaporation of the materials, we need to relocate this surface based on the criteria that evaporation takes place when temperature at any node reaches unity, and redistribute the meshes in the new domain. In order to establish the domain at every time step an iterative procedure, as illustrated in Fig. 4, is utilized. Figure 4 shows the iterative scheme to determine the interface location S. At any section a new position of interface location Z,,, = S is identified and calculated based on linear interpolation of two
M. J. KIMer al.
236
I
3
2
4
5 6 7 6 9 1011121314l51617
16
19 R
20
21 1
Fig. 3. Two-dimensional finite element mesh. other points between which the value of 8 crosses over. This position is expressed as below zpoin,= z3 + z,-s *4 _ & (8, - 1).
--X
--
____--___---
+-IY4
y*
Y4
1
Y
e,=i.i
y3
p
1.3
e,=
y2
f
1.4
e,=
Yl
____
(45)
Once Zpoinris identified, the length between 2 = S and ZE) is divided into equal lengths among the elements and then this new mesh distribution is used to obtain a new temperature distribution. This iteration procedure will be continued until a convergence
Y*
I
---p-
e, = 0.9
Ql
--
y3
-
_
* \
-
-I-
I
‘1
I
‘1
X
e4
93
6
j--+--t Y
Fig. 4. Iteration scheme.
I
1. t = 0.3 2. t = 0.4 3. t = 0.5 4. t = 0.6 5. 6. tt = = 0.7 0.8
-o.sl 4.0 -2.5 I
”
I
I
I
I
I
I
I
I
I
-1.0
-0.7
-0.3
0
0.3
0.7
1.0
2.5
4.0
Porition along the
length x
Fig. 5. The depth of the groove at different time t (LI = 0, N, = 0, Bi = 0.0001, Nk = 0.5).
PTP
237
Finite element modelling of the laser cutting process
-1.4 I -4.0
I -2.5
I -1.0
I I I I I -0.1 -0.3 0 0.3 0.7 Position along the length X
I 1.0
I 2.5
I 4.0
Fig. 6. The depth of the groove at different time I (CI = 0, N, = 0, Bi = 0.0001, Nk = 0.5).
in temperature at the surface is achieved within a certain assigned accuracy limit. RESULTS AND DISCUSSIONS
The finite element model is used to study the process of material removal from a solid surface subjected Gaussian laser beam of constant power. The evaporation rate, groove shape and temperature distribution in the media are studied for both steady state and transient conditions. Solutions are obtained for various operating parameters and cases with different boundary conditions, Transient state results Transient phenomena of material removal were studied first. Figures 5 and 6 show the variation of
depth of the groove along the axial position x with an increase in time for the case of U = 0. However, the rate of increment decreases as time increases and the heat losses are increased with an increase in area of interface. Figure 6 shows the final depth of the groove which remains steady with an increase in time. Figure 7 shows the variation of groove depth with an increase in time for the case with U = 1. The results show that the groove shape is similar to that in the case of U = 0, except it is slightly asymmetrical and reaches a stable position faster due to smaller residence time. Figures 8 and 9 show the comparison of groove shape given by both steady state and transient state analysis for the cases U = 0 and U = 1, respectively. It can be seen that this unsteady state solution is very close to the solution given by steady
m -0.6 t
-0.7-0.8 -“*9 -1.0 -1.1 -4.0
1.tm0.4 2. t = 0.5 3. t = 0.6 4. t = 0.7 5. t = 0.8 6. t = 4.5 I -2.5
I
I
-1.0
-0.7
6
,Vl -0.3
0
0.3
I 0.7
1.0
2.5
4.0
Position along the length x Fig. 7. The depth of the groove at different time r (U = I, N. = 0.01, Ei = 0.0001, Nk = 0.5). CAS 49,2-c
0 -0.1 -0.2 -0.3 -0.4 -0.5 -0.6 3
-0.7
*
-0.8 -0.9 -1.0 -1.1 -1.2
1.
unsteady state: t = 7.9
-1.3 -1.4 -4.0
-2.5
-1.0
-0.7
-0.3
0
0.3
0.7
1.0
2.5
10
Position along the length X
Fig. 8. The depth of the groove for steady state and unsteady state for U = 0, N, = 0, Bi = 0.0001, Nk = 0.5.
-0.1 -0.2 -0.3 -0.4 R
-0.5 -
:
-0.6 -0.7 -0.8 -0.9 -1.0 -1.1 -4.0
1. unsteady state: t =
4.5
2. steady state
I -2.5
I -1.0
I
I
I
I
I
I
-0.7
-0.3
0
0.3
0.7
1.0
Position along the
2.5
4.0
length X
Fig. 9. The depth of the groove for steady state and unsteady state for CJ= 1, N, = 0.01, Bi = 0.0001. Nk = 0.5.
B.C.:
B=Oatx=fq
Position along the length x Fig. 10. The depth of the groove for the different values of Nk. 238
Finite clement modelling of the laser cutting process
239
-0.8 -0.9 -
-1.0 _t.t _ -1.2 -
1. Bi 2. Bi 3. Bi 4. Bi
5 0.0001 re:0.01 = 0.1 = 0.3
5. Bi - 0.7 6. Bi = 0.9
-1.3 -4.0
-2.5
-1.0
-0.1
-0.3
0
0.3
0.7
1.0
2.5
4.0
Position along the length X Fig. 11. The depth of the groove for the different values of
state ana!ysis. The above results validate the usefulness of the steady state solutions. Therefore, the material removal rate and given shape can be predicted with good approximation by using steady state analysis.
Steady state results The steady state model is further used to study the material removal process with varying parameters and boundary conditions. Figures IO-16 show the results assuming the boundary condition at 2 = Z,, as 8 = 0. Figure 10 shows the variation of groove shapes with the different values of Nk. It can be seen from the definition of N,, that it is the ratio of conduction heat loss along the interface to the absorbed laser power. For large values of N,, the
Bi.
laser energy is mostly taken away by conductive heat losses and this causes the formation of the shallow groove. Decreasing Nk will cause a reduction in conductive heat losses and increases in the depth of the groove. However, as the groove becomes deeper, the radiation per unit surface area becomes smaller, resulting in less energy transfer into the material. At this stage, increasing Nk will not give higher material removal rate. In order to see the effects of surface convective heat losses on the groove shape and depth, results are presented in Fig. 11 for Nk = 0.5 with varying Biot numbers. For a fixed conduction resistance, a decrease in Biot number corresponds to increase in surface convection resistance, and this causes higher material removal rate. It can be seen from the results
B.C.=g=Oatr=fao,
-4.0
-2.3
-1.0
-0.7
-0.3
0
0.3
0.7
1.0
2.5
Porftion along the length x
Fig. 12. The depth of the groove for the different values of U.
1
-0.2 -0.4 ,+e+-+-+-+
-0.6 -0.8 -1.0 h -1.2 5 N -1.4 -1.6 -1.8 -2.0
x-x-x-.x-~
_w~xcx-x-x-x-x
-2.2
-1.0 -0.7
-0.3
0
0.3
0.7
1.0
2.5
4.0
Position along the length X Fig. 13. The distribution of the node points along the x-direction for U - 0, N, = 0, Nk = 0.5, Bi = 0.0001.
Position along the length X Fig. 14. The temperature distribution along the X direction for .!I = 0, N, = 0, Nk = 0.5, Bi = 0.0001. 0
-0.2 -0.4 -0.6 -0.8 -1.0 fi -1.2 5 N -1.4 -1.6 -1.8 -2.0
x-x-x-%-x_y~
_-x-x-x-x
-2.2
-0.7
-0.3
0
0.3
0.7
1.0
2.5
4.0
Position along the length x Fig. IS. The distribution of the node points along the X direction for U = 1, N, = 0.01, Nk = 0.5, Bi = 0.0001. 240
241
Finite element modelling of the laser cutting process
-1.0
-0.7
-0.3
0
0.3
0.7
1.0
1.5
4.0
Position along the length X 16. The temperature distribution along the X direction for U = 1, N, = 0.01, Nk = 0.5, Bi = 0.0001.
that groove depth and shape increase with a decrease in Biot number. Results also show that no appreciable change occurs in the groove if the Biot number is less than 0.01. Figure 12 shows the effect of nondimensional velocity U on the shape and material removal rate, shape of groove will change with the change of nondimensional velocity U, it relates the relative speed of the laser to that of thermal diffusion into the medium. For larger values of U material removal rate will be less because of small residence time. Therefore, both the groove depth and the slope of the left side edge will be decreased with an increase in value of U. Figure 13 shows the shape of groove and the positions of the finite element node points for the case with U = 0. As the material evaporates during the cutting process, meshes are regenerated in the media and nodes along each layer are displaced from their original position. Figure 14 shows the nodal temperature distribution along each layer. It can be seen that at the top layer the temperature of all nodes approach the evaporation temperature and this is uniform along throughout. At other layers, the nodes located away from the center of laser focal point have a lower temperature than those located at the center of the focal region. This is because less heat is taken away by conduction in the media than that used in evaporation for the material at the center node. Figures 15 and 16 show similar results for the case with U = 1. The groove shape and the temperature distribution is no longer symmetrical, and also less material is removed due to lower residence time. REFERENCES I. H. Soodak, Effects of heat transfer between gases and solids. Ph.D. thesis, Duke University, Durham, NC (1943).
2. H. G. Landau, Heat conduction in a melting solid. Q. appl Math. 8, 81-94 (1950). 3. J. E. Rogerson and G. A. Chayt, Total melting time in the ablating slab problem. J. appl. Phys. 42, 27 1I-27 13 (1971). 4. F. W. Dabby and U.-C. Paek, High-intensity laser-induced vaporization and explosion of solid material. IEEE J. Quantum Electronics QE-8, 106-l 11 (1972).
5. M. Von Allmen, Laser drilling velocity in metajs. J. appl. Phys. 47, 5460-5463
(1976).
6. H. E. Cline and T. R. Anthony, Heat treating and melting material with scanning laser or electron beam. J. appl. Phys. 48, 3895-3900 (1977). 7. J. Mazumder and W. M. Steen, Heat transfer model for CW laser material processing. J. appl. Phys. 51,941-947 (1980). 8. M. F. Modest and H. Abakians, Evaporative cutting of a semi-infinite body with a moving CW laser. Trans. ASME 108. 602-607 (1986). 9. S. C. Hsu, S. Chakravorty aid R: Mehrabian, Rapid melting and solidification of a surface layer. Metall. Trans.
9B, June (1978).
10. S. Biyikli and M. F. Modest, Effect of beam expansion
and focusing on evaporative cutting with a moving CW laser. J. Heal Transfer 110, 529-532 (1988). 11. C. Chan, J. Mazumdar and M. M. Chen, A two dimensional transient model for convection in laser melted pool. Mefall. Trans. lSA, 2175-2184 (1984). 12. C. P. Grigoropoulos, S. E. Long, A. F. Emery and W. E. D&her, Jr, Experimental and computational analysis of laser melting of thin silicon films. 1989 National Heat Transfer Conference, HTD-Vol. 113, Heat Transfer in Manufacturing and Materials Processing. 13. A. Minardi and P. J. Bishop, Temperature distribution within a metal subjected to irradiation by a laser of spatially varying intensity. 1989 National Heat Transfer Conference, HTD-Vol. 113, Heat Transfer in Manufacturing and Materials Processing. 14. J. M. Glass, H. P. Groger, R. J. Churchill, J. W. Lindau and T. E. Diller, Heat transfer in metallic glasses during laser cutting. 1989 National Heat Transfe; Conference, HTD-Vol. 113, Heat Transfer in Manufacturing and Materials Processing.