- Email: [email protected]
3D finite element analysis of evaporative laser cutting
Applied Mathematical Modelling 29 (2005) 938–954 www.elsevier.com/locate/apm
3D finite element analysis of evaporative laser cutting Meung Jung Kim
*
Department of Mechanical Engineering, Northern Illinois University, DeKalb, IL 60115, USA Received 1 June 2004; received in revised form 1 December 2004; accepted 8 February 2005 Available online 17 May 2005
Abstract A three-dimensional computational model of evaporative laser-cutting process has been developed using a finite element method. Steady heat transfer equation is used to model the laser-cutting process with a moving laser. The laser is assumed continuous wave Gaussian beam. The finite element surfaces on evaporation side are nonplanar and approximated by bilinear polynomial surfaces. Semi-infinite elements are introduced to approximate the semi-infinite domain. An iterative scheme is used to handle the geometric nonlinearity due to the unknown groove shape. The convergence studies are performed for various meshes. Numerical results about groove shapes and temperature distributions are presented and also compared with those by semi-analytical methods. 2005 Elsevier Inc. All rights reserved. Keywords: Finite element method; Evaporative laser cutting; Geometric nonlinearity; Groove shapes; Semi-infinite elements
1. Introduction The laser that was invented in 60s has found applications in many manufacturing processes primarily due to its precision process and high intensity [1–3]. The quality of the laser cut is of the utmost importance in laser processing because it would lead to an elimination of post-machining operations. Any improvement in laser cut quality would be of considerable significance. *
Tel.: +1 815 753 9965/9979; fax: +1 815 753 0416. E-mail address: [email protected]
0307-904X/$ - see front matter 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2005.02.015
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
Nomenclature Bi Biot number c specific heat heat of sublimation hig h convection heat transfer coefficient ^i; ^k unit vectors in the X- and Z-directions, respectively laser power density at the center of the beam I0 k thermal conductivity ^ n unit outward surface normal evaporation parameter Ne conduction parameter Nk conduction heat flux qk convection heat flux qh heat flux due to material evaporation qig heat flux due to laser radiation qL laser beam radius at the focal point Ro S(X, Y) groove depth s(x, y) non-dimensional groove depth final groove depth S1 T temperature ambient temperature T1 Tevap evaporation temperature t non-dimensional time u non-dimensional laser moving velocity U moving specimen or laser velocity x, y, z dimensionless spatial coordinates X, Y, Z spatial coordinates x1, x2, z1, z2 nodal coordinates of an element xF, yF, zF half the x, y, z-dimensions of the specimen xmin, xmax starting and ending x-coordinates of melting region on the specimen surface nodal z-coordinate at ith position zi new surface nodal z-coordinate at ith position znew i old ; z ! zactual actual nodal z-coordinate for iterative computation znew i i i Greek ao a e 1, g, n q h
letters absorptivity thermal diffusivity convergence limit for temperature and position dimensionless spatial coordinates of a field point density dimensionless temperature
939
940
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
hevap dimensionless evaporation temperature (hevap = 1) hi,m, hi,m+1 dimensionless temperatures at nodes zi,m and zi,m+1 k relaxation factor s time
High intensity laser beam can be directed to a narrow region in order to instantly evaporate material with very narrow heat affected zone. This ability to cut instantly with extremely narrow laser beam distinguishes it from other cutting methods. The key to success in precision cut by laser depends on many factors such as laser characteristics, material properties of the specimen, and manufacturing parameters. In precision manufacturing the quality of the cut is often measured based on the shape of the groove and amount of material removal. Therefore, better understanding of the process and thereby the roles of various parameters are essential to successful applications of laser-cutting process. There have been numerous investigations on laser applications. Some [4–8] investigated states of stresses in fracture, chemical compositions and properties, and heat transfer on different types of materials such as metals, composites, ceramics, and metallic glasses. Others [9–13] studied heat treatment effects on the material by laser irradiation. Considerable researches [14–18] have been done with heat transfer models on the effects of laser characteristics and material properties for the quality of laser processes. Other researches [1,19–21] are also found about pulsed lasers, melting of thin films, and reflections. As many applications need to take the melting into consideration, there are also many applications that rely on material evaporation such as cutting plastics and organic materials in medical operations. In addition, the current method for evaporative cutting can be extended to include melting pool in the future. This paper primarily focuses on the implementation of a three-dimensional finite element method for the first time in order to predict the groove shapes in evaporative laser cutting as an extension of the previous works by Kim et al. [22–25].
2. Mathematical formulation A typical laser cutting installation is shown in Fig. 1. The typical processes involved in evaaporative laser cutting are thermal in nature. When a laser beam strikes a material surface, several effects take place: reflection and absorption of the beam; conduction of heat into the material and loss of heat by convection and/or radiation from the material surface. The amount of energy absorbed and utilized in removing the material depends on the optical and thermo-physical properties of the material. The mathematical model describing the process of material removal from the surface subjected to high intensity laser beam can be found in Modest and Abakians [17,20] and Kim et al. [22–25]. It is assumed that there are three different regions on the surface subjected to laser beam as shown in Fig. 2. Region I is too far from the laser to have reached evaporation temperature, region II is the area in which evaporation takes place, and region III is the region in which evaporation has already taken place.
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
941
Fig. 1. Typical laser installation.
Fig. 2. Energy balance on the surface subject to laser.
Following assumptions are also made in deriving the model: (1) (2) (3) (4)
laser beam is of Gaussian type in a continuous mode, material moves at a constant relative velocity, material is isotropic and opaque with constant thermal and optical properties, material removal is a surface phenomenon and phase change from solid to vapor occurs in one step, (5) evaporated material is transparent and does not interfere with incident laser beam, (6) heat losses by convection and radiation from the surfaces to the environment can be approximated by using a single constant convection coefficient. Based on these assumptions, the mathematical statement of the problem can be written as follows: 2 oT o T o2 T o2 T ¼k þ þ ð1Þ qcU oX oX 2 oY 2 oZ 2 subjected to the boundary conditions at edges ^ qh qk ¼ ^
at X ¼ X F and Y ¼ Y F
and T ¼ T 0
at Z ¼ Z F for a finite model
ð2aÞ
942
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
or T ¼ T0
at X ¼ 1; Y ¼ 1
and Z ¼ 1 for a semi-infinite model
ð2bÞ
and the boundary condition on the surface subject to laser beam is obtained from the balance of heat transfer on the surface as ^ qk ¼ ^ qh þ ^ qig qL þ ^
ð3Þ ðX 2 þY 2 Þ=R2
where ^ qL ¼ ao I o ð^k ^ nÞe ;^ qh ¼ hðT T 1 Þ; ^qk ¼ kð^n rT Þ; ^qig ¼ qhig U ð^i ^nÞ; ^i and ^k are unit vectors in the X- and Z-direction, respectively, and ^n is the normal outward unit vector to surface. Here, qU ð^i ^ nÞ represents the rate of material removal when the specimen moves in negative X-direction with the speed U. With the introduction of the dimensionless variables as follows x¼
X ; Ro
y¼
qUhig Ne ¼ ; ao I o
Y ; Ro
z¼
Z ; Ro
sðx; yÞ ¼
kðT evap T 1 Þ Nk ¼ ; Ro ao I o
SðX ; Y Þ ; Ro
hRo ; Bi ¼ k
h¼
ðT T 1 Þ ; ðT evap T 1 Þ
URo u¼ ; a
ð4Þ
k a¼ qc
Eq. (1) and the boundary conditions (2) and (3) can be rewritten as: u
oh o2 h o2 h o2 h þ þ ¼ ox ox2 oy 2 oz2
ð5Þ
subjected to qk ¼ qh
at x ¼ xF and y ¼ y F
and h ¼ h0
h ¼ h1
at x ¼ 1; y ¼ 1 and z ¼ þ1
at z ¼ zF for a finite model
ð6aÞ
or ð6bÞ
and on the surface subject to laser Region I: z = 0, xF < x < xmin oh 2 2 N k Bih ¼ eðx þy Þ oz Region II: z = s(x, y), xmin < x < xmax " 2 2 #1=2 os oh os os ðx2 þy 2 Þ h ¼ 1; N e ¼ e N k Bih þ 1þ ox oz ox oy Region III: z = s1, xmax < x < xF " 2 #1=2 oh os1 2 2 N k Bih ¼ eðx þy Þ 1þ oz ox
ð7Þ
ð8Þ
ð9Þ
Here Bi is the Biot number representing the ratio of convection to conduction heat losses. u represents the ratio of relative speed of the work specimen to the thermal diffusivity of the material.
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
943
Ne is the ratio of energy utilized in evaporation of material and the absorbed laser energy while Nk represents the approximate ratio of conduction losses to the absorbed laser energy. The convection boundary condition in Eq. (2a) is now expressed in dimensionless form as qk ¼ Biðh h1 Þ
ð10Þ
The conduction term in Eq. (3) can now be expressed in terms of others in dimensionless form 2
qk ¼
ex ^ Ne ^ ðk ^ nÞ þ Biðh h1 Þ þ ði ^ nÞ Nk Nk
ð11Þ
Note that the regions I and III are subsets of region II and all regions can be handled by one type of region II in actual analysis.
3. Finite element formulation The variational formulation of the governing equation by weighted residual method leads to an integral form Z oh o2 h o2 h u 2 2 dh dX 0¼ ox ox oz X I Z oh oh odh oh odh oh oh u dh nx þ nz dX dh dC ¼ ox ox ox oz oz ox oz X C I Z oh oh odh oh odh u dh ð12Þ dX þ qk dh dC ¼ ox ox ox oz oz X C where X and C represent the domain and the boundary, respectively. The finite element formulation is obtained from this weak form of the variational formulation by introducing the shape functions /j(x, z) with nodal values hj on an elemental domain as h i ðeÞ ðeÞ K ij fhj g ¼ fF i g ð13Þ where ðeÞ K ij
ðeÞ
Fi
Z I o/j o/i o/j o/i o/j ¼ U /i Bi/i /j ds; þ þ dX ox ox ox oz oz Xe Ce I ¼ qk /i ds
i; j ¼ 1; 2; . . . ; N
Ce
Here Xe and Ce represent elemental domain and boundary. The dimensionless temperature h(x, z) is approximated by nodal values of temperature and shape functions as hðx; zÞ ¼
N X
hj /j ðx; zÞ
j¼1
and N is the number of node per element.
ð14Þ
944
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
By assembling the elemental finite element formulation the global finite element model of the problem can now be written as ½Kfhg ¼ fF g
ð15Þ
where [K] and {F} are known matrices and {h} is the unknown column vector. 4. Computational methods The governing equation is easier to express with moving specimen for the fixed frame to laser, but it is also easier to numerically implement with the moving laser for the fixed specimen. They are equivalent except the viewpoint and thus, in the present analysis the laser is handled as a moving source in positive X-direction with fixed specimen for easiness of numerical implementation below. Since the geometry (i.e., groove shape) is not known beforehand, computation begins with an assumed domain, which is the original shape of the specimen at the beginning. Once the temperature is calculated for the given domain, the nodal values of the surface temperatures in region II are examined if the boundary conditions in Eq. (8) are satisfied. If the temperature at any node is greater than the evaporation temperature, then the material at that node should have melted. In this paper a simple but effective linear interpolation as described in [24] and repeated below is used for the new position of a node. This simple scheme substantially reduces the computational time zi zbottom ¼ zi þ ðhi hevap Þ ð16Þ znew i hi hbottom is computed at node i, actual new value for next iteration is relaxed by Once znew i zactual ¼ ð1 kÞzold þ kznew i i i
ð17Þ
, is where k is a relaxation factor used to suppress oscillation in iteration. This new value, zactual i used to obtain a new domain. Since the nodes are moved in z-direction independently to simulate the material removal when the temperatures are greater than the melting temperature, the four surface nodes may not be planar and the surface integral due to the laser irradiation cannot be evaluated based on the planar assumption. Thus, the elemental surface with four nodes is approximated by a bilinear polynomial function as zðx; yÞ ¼ a0 þ a1 x þ a2 y þ a3 xy
ð18Þ
where the coefficients aiÕs can be obtained by imposing continuity conditions at four nodes as 8 9 2 38 9 a0 > z1 > 1 x1 y 1 x1 y 1 > > > > > > > > > > > > > 7 6 < z2 = 6 1 x2 y x2 y 7< a1 > = 2 2 7 ¼6 ð19Þ 6 1 x y x y 7> a > > z3 > > > 4 3 3 3 5> 2> 3 > > > > > > : > : > ; ; z4 1 x4 y 4 x4 y 4 a3 Here, ziÕs are nodal z-coordinates given by zi = z(xi, yi).
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
945
The region III where material has been removed is implemented by extending the maximum groove depths in the region II. Once the new shape of the domain is computed, the iteration continues until the relative sums of squared errors for both nodal temperatures and positions fall below a prescribed tolerance, e. " , #1=2 N N X X new old 2 new 2 ðhj hj Þ ðhj Þ ð20Þ Etemp ¼ j¼1
j¼1
"
Epos
N X 2 ¼ ðznew zold j j Þ
, #1=2 N X 2 ðznew j Þ
j¼1
ð21Þ
j¼1
Here N is the number of nodes in the domain. Even after the converged solution is obtained by iteration, it is possible that the temperature at a node has a value below the evaporation temperature but moved to a new position during iteration. If this happens, the node should be moved back to the original position and the iteration resumes. This causes numerical difficulty of unstable oscillation of errors during iteration.
5. Numerical results and discussion First, the three-dimensional results have been compared to two-dimensional results of previous works for various cases in Fig. 3. The dimensionless parameters for these cases with various speeds of laser are Bi ¼ 0.0001;
N e ¼ 0.001;
N k ¼ 0.4
ð22Þ
and the dimensions of the specimen are xF ¼ 8;
y F ¼ 8;
zF ¼ 2.5
ð23Þ
These values are chosen to roughly represent the specimen made of typical Aluminum cut by laser power of 1 kW with the beam focal radius of 0.1 mm that is subjected to natural convection by air. In this case the speed of laser beam in 1 m/s is converted to the dimensionless speed of 1 [25].
Fig. 3. Maximum groove depths of various cases in two- and three-dimensional cases.
946
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
The error criteria are e ¼ 0.01% for position error
ð24aÞ
e ¼ 0.1% for temperature error
ð24bÞ
The boundary conditions are given by Eq. (6a) and the laser is positioned at the center of the specimen. Fig. 3 shows the maximum groove depths against the laser velocity. Both two- and three-dimensional results converge with refined meshes, but the three-dimensional cases converge much faster than two-dimensional cases. The two-dimensional cases 40 · 5 and 80 · 5 or three-dimensional cases 40 · 40 · 5 and 80 · 80 · 5 meshes are very close and indistinguishable in the figure. It is noted that the values of the maximum groove depths of three-dimensional cases are approximately half of those in two-dimensional cases. This may be expected from the heat transfer to the third direction (positive and negative y directions) in three-dimensional cases that is absent in two-dimensional cases as illustrated in Fig. 4 (not to show actual direction of heat transfer that is normal to the surface). Also, the effect of number of Gauss integration points on maximum groove depths has been investigated for domain and surface integrals in (13). It was found (not shown here) that fine mesh of 40 · 40 · 5 and the minimum number of integration points, two, are good enough for numerically converged results. The second case considered deals with a semi-infinite body with the following parameters. U ¼ 1;
Bi ¼ 0.0001;
N e ¼ 0.01;
N k ¼ 0.1–0.005
ð25Þ
The semi-infinite elements used in this analysis are given in Appendix A. The boundary conditions are given by Eq. (6b). The effect of domain size and the mesh on maximum groove depths has been investigated. In computation, symmetry about x–z plane has been utilized in three-dimensional analyses. Table 1 shows two-dimensional results on the maximum groove depths as the mesh is refined and the domain is increased proportionally. The mesh of 20 · 15 elements in 2D with the domain size of 32 · 50 can be taken as converged results. Further, the result for domain sizes of 16 · 50 (not shown here) has been computed that is very close to the result for the case of 32 · 50. This suggests that the domain size Lx = 16 is good enough for accurate results.
x y
z
Fig. 4. Heat transfer in three directions at the cutting front.
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
947
Table 1 Effects of mesh and size of specimen on maximum groove depths with U = 1, Bi = 104, Ne = 102, Nk = 0.1 Lx · Lz
2D Nx · Nz
16 · 25
32 · 50
64 · 100
128 · 200
10 · 15
24.75272
20.55286
20.22619
20 · 15 40 · 15
14.55105 14.69527
21.52166 21.52356* 14.43881 14.43894
14.39967
Table 2 shows three-dimensional maximum groove depths as the mesh is doubled from left to right and the domain is increased proportionally to check the convergence with location of semiinfinite elements. Also, the number of elements in x- and y-directions from top to bottom shows the convergence. The mesh of 20 · 10 · 15 elements with the domain size of 32 · 16 · 50 yields reasonably converged results for this case. The values with * show very close results to preceding values indicating that domain size of Lx · Ly = 32 · 16 is good enough for semi-infinite dimensions in x- and y-directions. Modest et al. [17] predicted the max groove depth of roughly 10 in their semi-analytical analysis that is higher than the current numerical results. In Fig. 5 the max groove depths and temperature errors are shown during the iteration as well as the convergence in the error domain. In Fig. 6 the groove shape and temperature distribution with the half domain for the case 40 · 20 · 15 mesh in Table 2 are shown. Figs. 5(e) and 6(c) show that the temperature is more closely distributed in front of the moving laser showing the Doppler effect due to the moving source. It is interesting to see the caved-in contour plot at the bottom of the groove on the laser receding side. This suggests that the groove bottom cools down faster than the side surfaces. Further, the effect of number of integration points on maximum groove depths was also investigated for semi-infinite elements (not shown here). The results suggest that the specimen size of 16 and 8 in x- and y-directions is good for semi-infinite elements with two integration points. And consequently these values are used in the current analyses. Table 3 shows the case with Nk = 0.01. The convergence can be observed with 40 · 20 · 15 mesh. The present numerical results for maximum groove depth predict smaller value than the value 80 by Modest et al. [17]. Figs. 7 and 8 show typical changes of maximum groove depths, temperature and position errors, the temperature distribution, and the groove shape for the case of mesh 40 · 20 · 60 in Table 2 Effects of mesh and size of specimen on maximum groove depths with U = 1, Bi = 104, Ne = 102, Nk = 0.1 Lx · Ly · Lz
3D Starting Nx · Ny · Nz
*
16 · 8 · 25
32 · 16 · 50
64 · 32 · 100
10 · 5 · 15
7.61915
20 · 10 · 15 40 · 20 · 15
6.01683 5.80948
7.38836 7.38839* 5.765111 5.765111*
7.30045 7.30045* 5.74619*
The values are computed only increasing z-dimension from the preceding cases.
948
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
Fig. 5. (a) Maximum groove depth with iteration, (b) temperature error with iteration, (c) convergence behavior in error domain, (d) mid-plane groove shapes with iteration, and (e) mid-plane surface temperatures with iteration for the case of mesh 40 · 20 · 15 in Table 2.
Fig. 6. (a) 3D Groove shape and (b) 3D temperature distribution, (c) temperature contour plot on the groove surface (top view), and (d) heat flux in the mid-plane (at y = 0) for half domain for the case of mesh 40 · 20 · 15 in Table 2.
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
949
Table 3 Effects of mesh and size of specimen on maximum groove depths with U = 1, Bi = 104, Ne = 102, Nk = 0.01 Lx · Ly · Lz
3D Nx · Ny · Nz
*
16 · 8 · 200
32 · 16 · 400
64 · 32 · 800
10 · 5 · 15
69.19937
64.97447 64.97447*
20 · 10 · 15 40 · 20 · 15
47.52265 47.67044
66.09498 66.09512* 47.25097
The values in the table are computed only increasing z-dimension from the preceding cases.
Fig. 7. (a) Maximum groove depth with iteration, (b) temperature error with iteration, and (c) convergence behavior in error domain, (d) mid-plane groove shapes with iteration, and (e) mid-plane surface temperatures with iteration for the case of mesh 40 · 20 · 15 in Table 3.
Table 3. It is observed that the heat affected zone stretches far down from the laser with smaller conductivity of the material. In the following the sectional shapes of the groove are also presented. These shapes also compare well with those by Modest et al. [17] except the maximum groove depth. Here it is noted that Fig. 9(b) shows the backward groove depth for laser motion than forward groove for the moving specimen. A preliminary study with unsteady model (not shown here) shows forward groove with shallower depth in both cases of moving laser and moving specimen. Final case with Nk = 0.005 has been also studied and presented in Table 4 and Fig. 10. The converged results can be taken for the mesh 20 · 10 · 15 with domain size 32 · 16 · 400. The present
950
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
Fig. 8. (a) Groove shape and (b) temperature distribution, (c) temperature contour plot on the groove surface (top view), and (d) heat flux in the mid-plane (at y = 0) for half domain in the case of mesh 40 · 20 · 15 in Table 3.
Fig. 9. Groove section shapes (a) perpendicular to laser motion and (b) in the direction of laser motion for the case 40 · 20 · 15 in Table 3.
Table 4 Effects of mesh and size of specimen on maximum groove depths with U = 1, Bi = 104, Ne = 102, Nk = 0.005 Lx · Ly · Lz
3D Nx · Ny · Nz
*
16 · 8 · 300
32 · 16 · 400
64 · 32 · 800
10 · 5 · 15
125.77696
118.57616 118.57482*
20 · 0 · 15
74.76116
123.57431 123.57510* 74.90197 74.66565*
40 · 20 · 15
75.49031
The values are computed with same dimensions except doubled z-dimension.
numerical result for maximum groove depth predicts much smaller value that of 120 by Modest et al. [17] consistently. With the decrease of the conductivity most heat flux occurs along the laser motion horizontally.
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
951
Fig. 10. (a) Maximum groove depth with iteration, (b) temperature error with iteration, (c) error behavior in error domain, (d) groove shapes during iteration at mid-section (y = 0), (e) mid-section surface temperatures during iteration, (f) 3D groove shape, (g) 3D surface temperature, (h) temperature contour plot on the groove surface (top view), (i) midsection heat fluxes, (j) groove section shape in the direction of laser motion, and (k) groove section perpendicular to the laser motion for the case of mesh 40 · 20 · 15 in Table 4.
952
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
Fig. 10 shows similar behaviors of max groove depth with iteration as well as the groove shape and temperature distribution to the previous cases.
6. Conclusions A three-dimensional finite element model has been developed to analyze the evaporative lasercutting process based on steady heat conduction equation with constant laser velocity. The laser intensity is assumed to be sufficiently high to cause direct evaporation of the material from the surface of the medium. The laser side elemental surface is approximated by semi-quadratic polynomial function. Parametric study shows that the numerical results converge with the mesh refinement. The predicted groove shapes for a semi-infinite domain well compare with semi-analytical results by others except the max groove depths. The present analyses without any limiting assumptions predict smaller maximum groove depths than semi-analytical results. The temperature distributions show the heat-affected zone is not really limited close to the laser position. The Doppler effect is observed for a moving laser. The steady-state analyses show the backward groove shapes than forward ones in real laser-cutting process. The geometric nonlinearity due to the unknown groove shape has led to an iterative scheme that sometimes resulted in unstable oscillations during iterations.
Appendix A For eight node linear element (Fig. A.1), assuming that the element extends to infinity along zdirection the coordinates can be expressed as x¼
4 X j¼1
xj M j ;
y¼
4 X
yjM j;
z¼
j¼1
4 X
ðA:1Þ
zj M j
j¼1
where the mapping functions for coordinates are 1n 1g 2 ; 2 2 11 1þn 1þg 2 ; M3 ¼ 2 2 11 M1 ¼
1þn 1g 2 ; 2 2 11 1n 1þg 2 M4 ¼ 2 2 11 M2 ¼
ðA:2Þ
and the shape functions for field variable are standard linear shape functions. 4
3 η
8 5
1 ∞
7
2
ξ ζ
6
Fig. A.1. Natural coordinates of a semiinfinite element.
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
The Jacobian J in finite element formulation then becomes 3 3 2P 2 4 4 4 P P oM j oM j oM j ox oy oz x y z on j on j on j 7 6 on on on 7 6 j¼1 j¼1 j¼1 7 7 6 6 7 6 ox oy oz 7 6 4 4 4 P oM j P oM j P oM j 7 7 6 6 6 x y z 7 J ¼6 7¼ og j og j og j 7 6 og og og 7 6 j¼1 j¼1 j¼1 7 7 6 6 7 4 ox oy oz 5 6 4 4 4 5 4P P P oM j oM j oM j x y z o1 j o1 j o1 j o1 o1 o1 j¼1 j¼1 j¼1
953
ðA:3Þ
References [1] W.E. Lawson, Laser cutting of composites, in: Conference on Composites in Manufacturing, Technical Paper EM86-114, Los Angeles, 1986. [2] D.M. Roessler, Detroit looks to lasers, Mech. Eng. 112 (April 1990) (1986) 38–45. [3] J. Uhlenbusch, U. Bielesch, S. Klein, M. Napp, J.H. Schafer, Recent developments in metal processing with pulsed laser technology, Appl. Surf. Sci. 106 (1996) 228–234. [4] D.L. Carrol, J.A. Rothenflue, Experimental study of cutting thick aluminum and steel with a chemical oxygen– iodine laser using an N2 or O2 gas assist, J. Laser Appl. 9 (1997) 119–128. [5] S. Ghosh, B.P. Badgujar, G.L. Goswami, Parametric studies of cutting zircaloy-2 sheets with a laser beam, J. Laser Appl. 8 (1996) 143–148. [6] J.M. Glass et al., Heat transfer in metallic glasses during laser cutting, ASME Heat Transfer in Manufacturing and Materials Processing, New York, HTD Vol. 113, 1989, pp. 31–38. [7] J. Kusinski, Microstructure, chemical composition and properties of the surface layer of M2 steel after melting under different conditions, Appl. Surf. Sci. 86 (1995) 317–322. [8] K. Li, P. Sheng, Plane stress model for fracture of ceramics during laser cutting, Int. J. Mach. Tools Manufact. 35 (11) (1995) 1493–1506. [9] K. Brugger, Exact solutions for the temperature rise in a laser heated slab, J. Appl. Phys. 43 (1972) 577–583. [10] J.I. Masters, Problem of intense surface heating of a slab accompanied by change of phase, J. Appl. Phys. 27 (1956) 477–484. [11] A. Minardi, P.J. Bishop, Temperature distribution within a metal subjected to irradiation by a laser of spatially varying intensity, ASME Heat Transfer in Manufacturing and Materials Processing, New York, HTD Vol. 113, 1989, pp. 39–44. [12] N.N. Rykalin, A.A. Uglov, N.I. Makarov, Effects of peak frequency in a laser pulse on the heating of metal sheets, Sov. Phys.-Doklad. 12 (1967) 644–646. [13] R.E. Warren, M. Sparks, Laser heating of a slab having temperature-dependent surface absorptance, J. Appl. Phys. 50 (1979) 7952–7957. [14] B. Basu, J. Srinivasan, Numerical study of steady-state laser melting problem, Int. J. Heat Mass Transfer 31 (1988) 2331–2338. [15] K.A. Bunting, G. Cornfield, Toward a general theory of cutting: A relationship between the incident power density and the cut speed, ASME J. Heat Transfer 97 (1975) 116–121. [16] F.W. Dabby, U.C. Paek, High intensity laser-induced vaporization and explosion of solid material, IEEE J. Quant. Electron. QE-8 (1972) 106–111. [17] M.F. Modest, H. Abakian, Evaporative cutting of a semi-infinite body with a moving CW laser, ASME J. Heat Transfer 108 (1986) 597–601. [18] J.F. Reddy, Effects of High Power Laser Radiation, Academic Press, NY, 1971. [19] D. Maydan, Micromachining and image recording on thin films by laser beams, Bell Syst. Tech. J. 50 (1971) 1761–1789.
954
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
[20] M.F. Modest, H. Abakian, Heat conduction in a moving semi-infinite solid subjected to pulsed laser irradiation, ASME J. Heat Transfer 108 (1986) 602–607. [21] A.P. Zhuravel, A.G. Sivakov, O.G. Turutanov, I.M. Dmitrenko, A low temperature system with a pulse UV laser for scribing HTSC films and single crystals, Appl. Surf. Sci. 106 (1996) 321–325. [22] M.J. Kim, Z.H. Chen, P. Majumdar, Finite element modeling of the laser cutting process, J. Comput. Struct. 49 (1993) 231–241. [23] M.J. Kim, P. Majumdar, A computational model for high energy laser cutting process, Numer. Heat Transfer, Part A 27 (1995) 717–733. [24] M.J. Kim, Transient evaporative laser-cutting with boundary element method, J. Appl. Math. Modell. 25 (2000) 25–39. [25] M.J. Kim, Finite element analysis of evaporative cutting with a moving high energy pulsed laser, J. Appl. Math. Modell. 25 (2001) 203–220.
3D finite element analysis of evaporative laser cutting Meung Jung Kim
*
Department of Mechanical Engineering, Northern Illinois University, DeKalb, IL 60115, USA Received 1 June 2004; received in revised form 1 December 2004; accepted 8 February 2005 Available online 17 May 2005
Abstract A three-dimensional computational model of evaporative laser-cutting process has been developed using a finite element method. Steady heat transfer equation is used to model the laser-cutting process with a moving laser. The laser is assumed continuous wave Gaussian beam. The finite element surfaces on evaporation side are nonplanar and approximated by bilinear polynomial surfaces. Semi-infinite elements are introduced to approximate the semi-infinite domain. An iterative scheme is used to handle the geometric nonlinearity due to the unknown groove shape. The convergence studies are performed for various meshes. Numerical results about groove shapes and temperature distributions are presented and also compared with those by semi-analytical methods. 2005 Elsevier Inc. All rights reserved. Keywords: Finite element method; Evaporative laser cutting; Geometric nonlinearity; Groove shapes; Semi-infinite elements
1. Introduction The laser that was invented in 60s has found applications in many manufacturing processes primarily due to its precision process and high intensity [1–3]. The quality of the laser cut is of the utmost importance in laser processing because it would lead to an elimination of post-machining operations. Any improvement in laser cut quality would be of considerable significance. *
Tel.: +1 815 753 9965/9979; fax: +1 815 753 0416. E-mail address: [email protected]
0307-904X/$ - see front matter 2005 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2005.02.015
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
Nomenclature Bi Biot number c specific heat heat of sublimation hig h convection heat transfer coefficient ^i; ^k unit vectors in the X- and Z-directions, respectively laser power density at the center of the beam I0 k thermal conductivity ^ n unit outward surface normal evaporation parameter Ne conduction parameter Nk conduction heat flux qk convection heat flux qh heat flux due to material evaporation qig heat flux due to laser radiation qL laser beam radius at the focal point Ro S(X, Y) groove depth s(x, y) non-dimensional groove depth final groove depth S1 T temperature ambient temperature T1 Tevap evaporation temperature t non-dimensional time u non-dimensional laser moving velocity U moving specimen or laser velocity x, y, z dimensionless spatial coordinates X, Y, Z spatial coordinates x1, x2, z1, z2 nodal coordinates of an element xF, yF, zF half the x, y, z-dimensions of the specimen xmin, xmax starting and ending x-coordinates of melting region on the specimen surface nodal z-coordinate at ith position zi new surface nodal z-coordinate at ith position znew i old ; z ! zactual actual nodal z-coordinate for iterative computation znew i i i Greek ao a e 1, g, n q h
letters absorptivity thermal diffusivity convergence limit for temperature and position dimensionless spatial coordinates of a field point density dimensionless temperature
939
940
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
hevap dimensionless evaporation temperature (hevap = 1) hi,m, hi,m+1 dimensionless temperatures at nodes zi,m and zi,m+1 k relaxation factor s time
High intensity laser beam can be directed to a narrow region in order to instantly evaporate material with very narrow heat affected zone. This ability to cut instantly with extremely narrow laser beam distinguishes it from other cutting methods. The key to success in precision cut by laser depends on many factors such as laser characteristics, material properties of the specimen, and manufacturing parameters. In precision manufacturing the quality of the cut is often measured based on the shape of the groove and amount of material removal. Therefore, better understanding of the process and thereby the roles of various parameters are essential to successful applications of laser-cutting process. There have been numerous investigations on laser applications. Some [4–8] investigated states of stresses in fracture, chemical compositions and properties, and heat transfer on different types of materials such as metals, composites, ceramics, and metallic glasses. Others [9–13] studied heat treatment effects on the material by laser irradiation. Considerable researches [14–18] have been done with heat transfer models on the effects of laser characteristics and material properties for the quality of laser processes. Other researches [1,19–21] are also found about pulsed lasers, melting of thin films, and reflections. As many applications need to take the melting into consideration, there are also many applications that rely on material evaporation such as cutting plastics and organic materials in medical operations. In addition, the current method for evaporative cutting can be extended to include melting pool in the future. This paper primarily focuses on the implementation of a three-dimensional finite element method for the first time in order to predict the groove shapes in evaporative laser cutting as an extension of the previous works by Kim et al. [22–25].
2. Mathematical formulation A typical laser cutting installation is shown in Fig. 1. The typical processes involved in evaaporative laser cutting are thermal in nature. When a laser beam strikes a material surface, several effects take place: reflection and absorption of the beam; conduction of heat into the material and loss of heat by convection and/or radiation from the material surface. The amount of energy absorbed and utilized in removing the material depends on the optical and thermo-physical properties of the material. The mathematical model describing the process of material removal from the surface subjected to high intensity laser beam can be found in Modest and Abakians [17,20] and Kim et al. [22–25]. It is assumed that there are three different regions on the surface subjected to laser beam as shown in Fig. 2. Region I is too far from the laser to have reached evaporation temperature, region II is the area in which evaporation takes place, and region III is the region in which evaporation has already taken place.
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
941
Fig. 1. Typical laser installation.
Fig. 2. Energy balance on the surface subject to laser.
Following assumptions are also made in deriving the model: (1) (2) (3) (4)
laser beam is of Gaussian type in a continuous mode, material moves at a constant relative velocity, material is isotropic and opaque with constant thermal and optical properties, material removal is a surface phenomenon and phase change from solid to vapor occurs in one step, (5) evaporated material is transparent and does not interfere with incident laser beam, (6) heat losses by convection and radiation from the surfaces to the environment can be approximated by using a single constant convection coefficient. Based on these assumptions, the mathematical statement of the problem can be written as follows: 2 oT o T o2 T o2 T ¼k þ þ ð1Þ qcU oX oX 2 oY 2 oZ 2 subjected to the boundary conditions at edges ^ qh qk ¼ ^
at X ¼ X F and Y ¼ Y F
and T ¼ T 0
at Z ¼ Z F for a finite model
ð2aÞ
942
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
or T ¼ T0
at X ¼ 1; Y ¼ 1
and Z ¼ 1 for a semi-infinite model
ð2bÞ
and the boundary condition on the surface subject to laser beam is obtained from the balance of heat transfer on the surface as ^ qk ¼ ^ qh þ ^ qig qL þ ^
ð3Þ ðX 2 þY 2 Þ=R2
where ^ qL ¼ ao I o ð^k ^ nÞe ;^ qh ¼ hðT T 1 Þ; ^qk ¼ kð^n rT Þ; ^qig ¼ qhig U ð^i ^nÞ; ^i and ^k are unit vectors in the X- and Z-direction, respectively, and ^n is the normal outward unit vector to surface. Here, qU ð^i ^ nÞ represents the rate of material removal when the specimen moves in negative X-direction with the speed U. With the introduction of the dimensionless variables as follows x¼
X ; Ro
y¼
qUhig Ne ¼ ; ao I o
Y ; Ro
z¼
Z ; Ro
sðx; yÞ ¼
kðT evap T 1 Þ Nk ¼ ; Ro ao I o
SðX ; Y Þ ; Ro
hRo ; Bi ¼ k
h¼
ðT T 1 Þ ; ðT evap T 1 Þ
URo u¼ ; a
ð4Þ
k a¼ qc
Eq. (1) and the boundary conditions (2) and (3) can be rewritten as: u
oh o2 h o2 h o2 h þ þ ¼ ox ox2 oy 2 oz2
ð5Þ
subjected to qk ¼ qh
at x ¼ xF and y ¼ y F
and h ¼ h0
h ¼ h1
at x ¼ 1; y ¼ 1 and z ¼ þ1
at z ¼ zF for a finite model
ð6aÞ
or ð6bÞ
and on the surface subject to laser Region I: z = 0, xF < x < xmin oh 2 2 N k Bih ¼ eðx þy Þ oz Region II: z = s(x, y), xmin < x < xmax " 2 2 #1=2 os oh os os ðx2 þy 2 Þ h ¼ 1; N e ¼ e N k Bih þ 1þ ox oz ox oy Region III: z = s1, xmax < x < xF " 2 #1=2 oh os1 2 2 N k Bih ¼ eðx þy Þ 1þ oz ox
ð7Þ
ð8Þ
ð9Þ
Here Bi is the Biot number representing the ratio of convection to conduction heat losses. u represents the ratio of relative speed of the work specimen to the thermal diffusivity of the material.
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
943
Ne is the ratio of energy utilized in evaporation of material and the absorbed laser energy while Nk represents the approximate ratio of conduction losses to the absorbed laser energy. The convection boundary condition in Eq. (2a) is now expressed in dimensionless form as qk ¼ Biðh h1 Þ
ð10Þ
The conduction term in Eq. (3) can now be expressed in terms of others in dimensionless form 2
qk ¼
ex ^ Ne ^ ðk ^ nÞ þ Biðh h1 Þ þ ði ^ nÞ Nk Nk
ð11Þ
Note that the regions I and III are subsets of region II and all regions can be handled by one type of region II in actual analysis.
3. Finite element formulation The variational formulation of the governing equation by weighted residual method leads to an integral form Z oh o2 h o2 h u 2 2 dh dX 0¼ ox ox oz X I Z oh oh odh oh odh oh oh u dh nx þ nz dX dh dC ¼ ox ox ox oz oz ox oz X C I Z oh oh odh oh odh u dh ð12Þ dX þ qk dh dC ¼ ox ox ox oz oz X C where X and C represent the domain and the boundary, respectively. The finite element formulation is obtained from this weak form of the variational formulation by introducing the shape functions /j(x, z) with nodal values hj on an elemental domain as h i ðeÞ ðeÞ K ij fhj g ¼ fF i g ð13Þ where ðeÞ K ij
ðeÞ
Fi
Z I o/j o/i o/j o/i o/j ¼ U /i Bi/i /j ds; þ þ dX ox ox ox oz oz Xe Ce I ¼ qk /i ds
i; j ¼ 1; 2; . . . ; N
Ce
Here Xe and Ce represent elemental domain and boundary. The dimensionless temperature h(x, z) is approximated by nodal values of temperature and shape functions as hðx; zÞ ¼
N X
hj /j ðx; zÞ
j¼1
and N is the number of node per element.
ð14Þ
944
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
By assembling the elemental finite element formulation the global finite element model of the problem can now be written as ½Kfhg ¼ fF g
ð15Þ
where [K] and {F} are known matrices and {h} is the unknown column vector. 4. Computational methods The governing equation is easier to express with moving specimen for the fixed frame to laser, but it is also easier to numerically implement with the moving laser for the fixed specimen. They are equivalent except the viewpoint and thus, in the present analysis the laser is handled as a moving source in positive X-direction with fixed specimen for easiness of numerical implementation below. Since the geometry (i.e., groove shape) is not known beforehand, computation begins with an assumed domain, which is the original shape of the specimen at the beginning. Once the temperature is calculated for the given domain, the nodal values of the surface temperatures in region II are examined if the boundary conditions in Eq. (8) are satisfied. If the temperature at any node is greater than the evaporation temperature, then the material at that node should have melted. In this paper a simple but effective linear interpolation as described in [24] and repeated below is used for the new position of a node. This simple scheme substantially reduces the computational time zi zbottom ¼ zi þ ðhi hevap Þ ð16Þ znew i hi hbottom is computed at node i, actual new value for next iteration is relaxed by Once znew i zactual ¼ ð1 kÞzold þ kznew i i i
ð17Þ
, is where k is a relaxation factor used to suppress oscillation in iteration. This new value, zactual i used to obtain a new domain. Since the nodes are moved in z-direction independently to simulate the material removal when the temperatures are greater than the melting temperature, the four surface nodes may not be planar and the surface integral due to the laser irradiation cannot be evaluated based on the planar assumption. Thus, the elemental surface with four nodes is approximated by a bilinear polynomial function as zðx; yÞ ¼ a0 þ a1 x þ a2 y þ a3 xy
ð18Þ
where the coefficients aiÕs can be obtained by imposing continuity conditions at four nodes as 8 9 2 38 9 a0 > z1 > 1 x1 y 1 x1 y 1 > > > > > > > > > > > > > 7 6 < z2 = 6 1 x2 y x2 y 7< a1 > = 2 2 7 ¼6 ð19Þ 6 1 x y x y 7> a > > z3 > > > 4 3 3 3 5> 2> 3 > > > > > > : > : > ; ; z4 1 x4 y 4 x4 y 4 a3 Here, ziÕs are nodal z-coordinates given by zi = z(xi, yi).
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
945
The region III where material has been removed is implemented by extending the maximum groove depths in the region II. Once the new shape of the domain is computed, the iteration continues until the relative sums of squared errors for both nodal temperatures and positions fall below a prescribed tolerance, e. " , #1=2 N N X X new old 2 new 2 ðhj hj Þ ðhj Þ ð20Þ Etemp ¼ j¼1
j¼1
"
Epos
N X 2 ¼ ðznew zold j j Þ
, #1=2 N X 2 ðznew j Þ
j¼1
ð21Þ
j¼1
Here N is the number of nodes in the domain. Even after the converged solution is obtained by iteration, it is possible that the temperature at a node has a value below the evaporation temperature but moved to a new position during iteration. If this happens, the node should be moved back to the original position and the iteration resumes. This causes numerical difficulty of unstable oscillation of errors during iteration.
5. Numerical results and discussion First, the three-dimensional results have been compared to two-dimensional results of previous works for various cases in Fig. 3. The dimensionless parameters for these cases with various speeds of laser are Bi ¼ 0.0001;
N e ¼ 0.001;
N k ¼ 0.4
ð22Þ
and the dimensions of the specimen are xF ¼ 8;
y F ¼ 8;
zF ¼ 2.5
ð23Þ
These values are chosen to roughly represent the specimen made of typical Aluminum cut by laser power of 1 kW with the beam focal radius of 0.1 mm that is subjected to natural convection by air. In this case the speed of laser beam in 1 m/s is converted to the dimensionless speed of 1 [25].
Fig. 3. Maximum groove depths of various cases in two- and three-dimensional cases.
946
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
The error criteria are e ¼ 0.01% for position error
ð24aÞ
e ¼ 0.1% for temperature error
ð24bÞ
The boundary conditions are given by Eq. (6a) and the laser is positioned at the center of the specimen. Fig. 3 shows the maximum groove depths against the laser velocity. Both two- and three-dimensional results converge with refined meshes, but the three-dimensional cases converge much faster than two-dimensional cases. The two-dimensional cases 40 · 5 and 80 · 5 or three-dimensional cases 40 · 40 · 5 and 80 · 80 · 5 meshes are very close and indistinguishable in the figure. It is noted that the values of the maximum groove depths of three-dimensional cases are approximately half of those in two-dimensional cases. This may be expected from the heat transfer to the third direction (positive and negative y directions) in three-dimensional cases that is absent in two-dimensional cases as illustrated in Fig. 4 (not to show actual direction of heat transfer that is normal to the surface). Also, the effect of number of Gauss integration points on maximum groove depths has been investigated for domain and surface integrals in (13). It was found (not shown here) that fine mesh of 40 · 40 · 5 and the minimum number of integration points, two, are good enough for numerically converged results. The second case considered deals with a semi-infinite body with the following parameters. U ¼ 1;
Bi ¼ 0.0001;
N e ¼ 0.01;
N k ¼ 0.1–0.005
ð25Þ
The semi-infinite elements used in this analysis are given in Appendix A. The boundary conditions are given by Eq. (6b). The effect of domain size and the mesh on maximum groove depths has been investigated. In computation, symmetry about x–z plane has been utilized in three-dimensional analyses. Table 1 shows two-dimensional results on the maximum groove depths as the mesh is refined and the domain is increased proportionally. The mesh of 20 · 15 elements in 2D with the domain size of 32 · 50 can be taken as converged results. Further, the result for domain sizes of 16 · 50 (not shown here) has been computed that is very close to the result for the case of 32 · 50. This suggests that the domain size Lx = 16 is good enough for accurate results.
x y
z
Fig. 4. Heat transfer in three directions at the cutting front.
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
947
Table 1 Effects of mesh and size of specimen on maximum groove depths with U = 1, Bi = 104, Ne = 102, Nk = 0.1 Lx · Lz
2D Nx · Nz
16 · 25
32 · 50
64 · 100
128 · 200
10 · 15
24.75272
20.55286
20.22619
20 · 15 40 · 15
14.55105 14.69527
21.52166 21.52356* 14.43881 14.43894
14.39967
Table 2 shows three-dimensional maximum groove depths as the mesh is doubled from left to right and the domain is increased proportionally to check the convergence with location of semiinfinite elements. Also, the number of elements in x- and y-directions from top to bottom shows the convergence. The mesh of 20 · 10 · 15 elements with the domain size of 32 · 16 · 50 yields reasonably converged results for this case. The values with * show very close results to preceding values indicating that domain size of Lx · Ly = 32 · 16 is good enough for semi-infinite dimensions in x- and y-directions. Modest et al. [17] predicted the max groove depth of roughly 10 in their semi-analytical analysis that is higher than the current numerical results. In Fig. 5 the max groove depths and temperature errors are shown during the iteration as well as the convergence in the error domain. In Fig. 6 the groove shape and temperature distribution with the half domain for the case 40 · 20 · 15 mesh in Table 2 are shown. Figs. 5(e) and 6(c) show that the temperature is more closely distributed in front of the moving laser showing the Doppler effect due to the moving source. It is interesting to see the caved-in contour plot at the bottom of the groove on the laser receding side. This suggests that the groove bottom cools down faster than the side surfaces. Further, the effect of number of integration points on maximum groove depths was also investigated for semi-infinite elements (not shown here). The results suggest that the specimen size of 16 and 8 in x- and y-directions is good for semi-infinite elements with two integration points. And consequently these values are used in the current analyses. Table 3 shows the case with Nk = 0.01. The convergence can be observed with 40 · 20 · 15 mesh. The present numerical results for maximum groove depth predict smaller value than the value 80 by Modest et al. [17]. Figs. 7 and 8 show typical changes of maximum groove depths, temperature and position errors, the temperature distribution, and the groove shape for the case of mesh 40 · 20 · 60 in Table 2 Effects of mesh and size of specimen on maximum groove depths with U = 1, Bi = 104, Ne = 102, Nk = 0.1 Lx · Ly · Lz
3D Starting Nx · Ny · Nz
*
16 · 8 · 25
32 · 16 · 50
64 · 32 · 100
10 · 5 · 15
7.61915
20 · 10 · 15 40 · 20 · 15
6.01683 5.80948
7.38836 7.38839* 5.765111 5.765111*
7.30045 7.30045* 5.74619*
The values are computed only increasing z-dimension from the preceding cases.
948
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
Fig. 5. (a) Maximum groove depth with iteration, (b) temperature error with iteration, (c) convergence behavior in error domain, (d) mid-plane groove shapes with iteration, and (e) mid-plane surface temperatures with iteration for the case of mesh 40 · 20 · 15 in Table 2.
Fig. 6. (a) 3D Groove shape and (b) 3D temperature distribution, (c) temperature contour plot on the groove surface (top view), and (d) heat flux in the mid-plane (at y = 0) for half domain for the case of mesh 40 · 20 · 15 in Table 2.
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
949
Table 3 Effects of mesh and size of specimen on maximum groove depths with U = 1, Bi = 104, Ne = 102, Nk = 0.01 Lx · Ly · Lz
3D Nx · Ny · Nz
*
16 · 8 · 200
32 · 16 · 400
64 · 32 · 800
10 · 5 · 15
69.19937
64.97447 64.97447*
20 · 10 · 15 40 · 20 · 15
47.52265 47.67044
66.09498 66.09512* 47.25097
The values in the table are computed only increasing z-dimension from the preceding cases.
Fig. 7. (a) Maximum groove depth with iteration, (b) temperature error with iteration, and (c) convergence behavior in error domain, (d) mid-plane groove shapes with iteration, and (e) mid-plane surface temperatures with iteration for the case of mesh 40 · 20 · 15 in Table 3.
Table 3. It is observed that the heat affected zone stretches far down from the laser with smaller conductivity of the material. In the following the sectional shapes of the groove are also presented. These shapes also compare well with those by Modest et al. [17] except the maximum groove depth. Here it is noted that Fig. 9(b) shows the backward groove depth for laser motion than forward groove for the moving specimen. A preliminary study with unsteady model (not shown here) shows forward groove with shallower depth in both cases of moving laser and moving specimen. Final case with Nk = 0.005 has been also studied and presented in Table 4 and Fig. 10. The converged results can be taken for the mesh 20 · 10 · 15 with domain size 32 · 16 · 400. The present
950
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
Fig. 8. (a) Groove shape and (b) temperature distribution, (c) temperature contour plot on the groove surface (top view), and (d) heat flux in the mid-plane (at y = 0) for half domain in the case of mesh 40 · 20 · 15 in Table 3.
Fig. 9. Groove section shapes (a) perpendicular to laser motion and (b) in the direction of laser motion for the case 40 · 20 · 15 in Table 3.
Table 4 Effects of mesh and size of specimen on maximum groove depths with U = 1, Bi = 104, Ne = 102, Nk = 0.005 Lx · Ly · Lz
3D Nx · Ny · Nz
*
16 · 8 · 300
32 · 16 · 400
64 · 32 · 800
10 · 5 · 15
125.77696
118.57616 118.57482*
20 · 0 · 15
74.76116
123.57431 123.57510* 74.90197 74.66565*
40 · 20 · 15
75.49031
The values are computed with same dimensions except doubled z-dimension.
numerical result for maximum groove depth predicts much smaller value that of 120 by Modest et al. [17] consistently. With the decrease of the conductivity most heat flux occurs along the laser motion horizontally.
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
951
Fig. 10. (a) Maximum groove depth with iteration, (b) temperature error with iteration, (c) error behavior in error domain, (d) groove shapes during iteration at mid-section (y = 0), (e) mid-section surface temperatures during iteration, (f) 3D groove shape, (g) 3D surface temperature, (h) temperature contour plot on the groove surface (top view), (i) midsection heat fluxes, (j) groove section shape in the direction of laser motion, and (k) groove section perpendicular to the laser motion for the case of mesh 40 · 20 · 15 in Table 4.
952
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
Fig. 10 shows similar behaviors of max groove depth with iteration as well as the groove shape and temperature distribution to the previous cases.
6. Conclusions A three-dimensional finite element model has been developed to analyze the evaporative lasercutting process based on steady heat conduction equation with constant laser velocity. The laser intensity is assumed to be sufficiently high to cause direct evaporation of the material from the surface of the medium. The laser side elemental surface is approximated by semi-quadratic polynomial function. Parametric study shows that the numerical results converge with the mesh refinement. The predicted groove shapes for a semi-infinite domain well compare with semi-analytical results by others except the max groove depths. The present analyses without any limiting assumptions predict smaller maximum groove depths than semi-analytical results. The temperature distributions show the heat-affected zone is not really limited close to the laser position. The Doppler effect is observed for a moving laser. The steady-state analyses show the backward groove shapes than forward ones in real laser-cutting process. The geometric nonlinearity due to the unknown groove shape has led to an iterative scheme that sometimes resulted in unstable oscillations during iterations.
Appendix A For eight node linear element (Fig. A.1), assuming that the element extends to infinity along zdirection the coordinates can be expressed as x¼
4 X j¼1
xj M j ;
y¼
4 X
yjM j;
z¼
j¼1
4 X
ðA:1Þ
zj M j
j¼1
where the mapping functions for coordinates are 1n 1g 2 ; 2 2 11 1þn 1þg 2 ; M3 ¼ 2 2 11 M1 ¼
1þn 1g 2 ; 2 2 11 1n 1þg 2 M4 ¼ 2 2 11 M2 ¼
ðA:2Þ
and the shape functions for field variable are standard linear shape functions. 4
3 η
8 5
1 ∞
7
2
ξ ζ
6
Fig. A.1. Natural coordinates of a semiinfinite element.
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
The Jacobian J in finite element formulation then becomes 3 3 2P 2 4 4 4 P P oM j oM j oM j ox oy oz x y z on j on j on j 7 6 on on on 7 6 j¼1 j¼1 j¼1 7 7 6 6 7 6 ox oy oz 7 6 4 4 4 P oM j P oM j P oM j 7 7 6 6 6 x y z 7 J ¼6 7¼ og j og j og j 7 6 og og og 7 6 j¼1 j¼1 j¼1 7 7 6 6 7 4 ox oy oz 5 6 4 4 4 5 4P P P oM j oM j oM j x y z o1 j o1 j o1 j o1 o1 o1 j¼1 j¼1 j¼1
953
ðA:3Þ
References [1] W.E. Lawson, Laser cutting of composites, in: Conference on Composites in Manufacturing, Technical Paper EM86-114, Los Angeles, 1986. [2] D.M. Roessler, Detroit looks to lasers, Mech. Eng. 112 (April 1990) (1986) 38–45. [3] J. Uhlenbusch, U. Bielesch, S. Klein, M. Napp, J.H. Schafer, Recent developments in metal processing with pulsed laser technology, Appl. Surf. Sci. 106 (1996) 228–234. [4] D.L. Carrol, J.A. Rothenflue, Experimental study of cutting thick aluminum and steel with a chemical oxygen– iodine laser using an N2 or O2 gas assist, J. Laser Appl. 9 (1997) 119–128. [5] S. Ghosh, B.P. Badgujar, G.L. Goswami, Parametric studies of cutting zircaloy-2 sheets with a laser beam, J. Laser Appl. 8 (1996) 143–148. [6] J.M. Glass et al., Heat transfer in metallic glasses during laser cutting, ASME Heat Transfer in Manufacturing and Materials Processing, New York, HTD Vol. 113, 1989, pp. 31–38. [7] J. Kusinski, Microstructure, chemical composition and properties of the surface layer of M2 steel after melting under different conditions, Appl. Surf. Sci. 86 (1995) 317–322. [8] K. Li, P. Sheng, Plane stress model for fracture of ceramics during laser cutting, Int. J. Mach. Tools Manufact. 35 (11) (1995) 1493–1506. [9] K. Brugger, Exact solutions for the temperature rise in a laser heated slab, J. Appl. Phys. 43 (1972) 577–583. [10] J.I. Masters, Problem of intense surface heating of a slab accompanied by change of phase, J. Appl. Phys. 27 (1956) 477–484. [11] A. Minardi, P.J. Bishop, Temperature distribution within a metal subjected to irradiation by a laser of spatially varying intensity, ASME Heat Transfer in Manufacturing and Materials Processing, New York, HTD Vol. 113, 1989, pp. 39–44. [12] N.N. Rykalin, A.A. Uglov, N.I. Makarov, Effects of peak frequency in a laser pulse on the heating of metal sheets, Sov. Phys.-Doklad. 12 (1967) 644–646. [13] R.E. Warren, M. Sparks, Laser heating of a slab having temperature-dependent surface absorptance, J. Appl. Phys. 50 (1979) 7952–7957. [14] B. Basu, J. Srinivasan, Numerical study of steady-state laser melting problem, Int. J. Heat Mass Transfer 31 (1988) 2331–2338. [15] K.A. Bunting, G. Cornfield, Toward a general theory of cutting: A relationship between the incident power density and the cut speed, ASME J. Heat Transfer 97 (1975) 116–121. [16] F.W. Dabby, U.C. Paek, High intensity laser-induced vaporization and explosion of solid material, IEEE J. Quant. Electron. QE-8 (1972) 106–111. [17] M.F. Modest, H. Abakian, Evaporative cutting of a semi-infinite body with a moving CW laser, ASME J. Heat Transfer 108 (1986) 597–601. [18] J.F. Reddy, Effects of High Power Laser Radiation, Academic Press, NY, 1971. [19] D. Maydan, Micromachining and image recording on thin films by laser beams, Bell Syst. Tech. J. 50 (1971) 1761–1789.
954
M.J. Kim / Applied Mathematical Modelling 29 (2005) 938–954
[20] M.F. Modest, H. Abakian, Heat conduction in a moving semi-infinite solid subjected to pulsed laser irradiation, ASME J. Heat Transfer 108 (1986) 602–607. [21] A.P. Zhuravel, A.G. Sivakov, O.G. Turutanov, I.M. Dmitrenko, A low temperature system with a pulse UV laser for scribing HTSC films and single crystals, Appl. Surf. Sci. 106 (1996) 321–325. [22] M.J. Kim, Z.H. Chen, P. Majumdar, Finite element modeling of the laser cutting process, J. Comput. Struct. 49 (1993) 231–241. [23] M.J. Kim, P. Majumdar, A computational model for high energy laser cutting process, Numer. Heat Transfer, Part A 27 (1995) 717–733. [24] M.J. Kim, Transient evaporative laser-cutting with boundary element method, J. Appl. Math. Modell. 25 (2000) 25–39. [25] M.J. Kim, Finite element analysis of evaporative cutting with a moving high energy pulsed laser, J. Appl. Math. Modell. 25 (2001) 203–220.