Effect of Process Parameters and Material Properties on Laser Micromachining of Microchannels

Laser micromachining has emerged as a promising technique for mass production of microfluidic devices. However, control and optimization of process parameters, and design of substrate materials are still ongoing challenges for the widespread application of laser micromachining. This article reports a systematic study on the effect of laser system parameters and thermo-physical properties of substrate materials on laser micromachining. Three dimensional transient heat conduction equation with a Gaussian laser heat source was solved using finite element based Multiphysics software COMSOL 5.2a. Large heat convection coefficients were used to consider the rapid phase transition of the material during the laser treatment. The depth of the laser cut was measured by removing material at a pre-set temperature. The grid independent analysis was performed for ensuring the accuracy of the model. The results show that laser power and scanning speed have a strong effect on the channel depth, while the level of focus of the laser beam contributes in determining both the depth and width of the channel. Higher thermal conductivity results deeper in cuts, in contrast the higher specific heat produces shallower channels for a given condition. These findings can help in designing and optimizing process parameters for laser micromachining of microfluidic devices.


Introduction
Microfluidic devices have been shown to have a broad range of potential applications in chemistry, biology, engineering, and biomedical sciences such as polymerase chain reaction (PCR) [1][2][3], protein separation and analysis [4][5][6], sample mixing [7], diagnostics [8], and chemical reactions [7,9] etc. The use of microfluidics offers several competitive advantages such as requiring smaller sample volumes, fewer reagents, greater portability, faster processing, easier disposability, and better compatibility with smart devices [8]. However, mass production of these devices without compromise on performance is an ongoing challenge.
Photolithography, chemical etching, soft lithography, hot embossing and injection molding are the most commonly used microfabrication processes for fabricating microfluidic devices [10]. Many of these processes are either expensive, incompatible for mass production, or limited to use on certain materials as substrates [11][12][13]. For instance, chemical etching and photolithographic process mainly use silica based materials such as silicon wafer, glass and polydimethylsiloxane (PDMS) as substrates [2]. Recently laser (Burlington, MA, USA). The transient heat conduction equation was solved considering a continuous, Gaussian laser beam as a moving heat source. The effect of laser spot size, laser scanning speed and laser power on the resulting channel width, depth, and shape were investigated. A parametric study of the effect of various relevant thermos-physical properties of the substrate materials was also investigated. The rest of the paper is organized as follows: Section 2 contains theory, model formulation, boundary conditions, convergence studies; Section 3 provides results and discussion; and finally, Section 4 provides conclusions of the study.

Theory
When a continuous IR laser beam, such as the CO 2 laser, is moved across the surface of a polymer, the laser energy is absorbed as heat, which rapidly increases the temperature of the material in the irradiated region. At the spot where the laser meets the surface, a portion of the material is vaporized, while another portion melts into a liquid form. The vaporized material expands and escapes, driving the liquid polymer outward along with it. In this way, some of the liquid material is ejected, while the rest remains in the channel and re-solidifies in the wake of the beam. By moving the laser across the surface, a channel and other similar structures can be cut into the material. The cross section of the resulting channel depends on the material which is mostly regulated by thermal diffusivity of that material, as well as various laser parameters. Therefore, the CO 2 laser micromachining of microfluidic channels can be modeled as a thermal event [33] using the transient heat equation.
Although the material undergoes a phase change as it is heated by the laser, most of the heating process occurs in the solid phase. Once the material reaches the point of phase change, it quickly changes phase and is removed from the bulk solid material. Therefore, the phase change and subsequent material removal can be modeled as a convection boundary condition in a transient heat conduction equation with laser irradiation as a moving heat source.

Governing Equations
To model laser cutting as a thermal event, the transient heat equation is given by [34]: where ρ is the density, C p is the constant pressure specific heat, and k is the thermal conductivity of the material. The . q V term represents volumetric heat generation due to the absorption of the laser beam at the surface. The specific laser energy absorption is given by using Beer's Law [35]: where α is the material's absorption coefficient (m −1 ) and I(r,z) is the intensity of the laser at a distance r from the axis of the beam, the z-coordinate is the distance from the surface of the substrate measured downwards in the direction of the beam. The intensity profile of the laser beams is given by a Gaussian function [36]: where I(0,z) is the peak intensity along the beam axis, I(0, z) = 2P 0 πw(z+d) 2 , P 0 is the total power of the laser beam, w(z) is the beam radius, or spot size, at which point the beam intensity falls to 1/e 2 of its value along the beam axis, and d is the distance from the substrate to the focal point of the focusing lens as shown in Figure 1a. The size of the beam, described by the beam radius, is given by the following function: where λ is the wavelength of the laser and w 0 is the waist radius, the minimum radius of the laser beam. Assuming the beam radius is approximately constant over relevant z values, the beam radius becomes a constant w 0 . Combining and simplifying Equations (2)-(4), and using r 2 = x 2 + y 2 , For laser cutting, the laser beam is moved continuously across the surface. Assuming the laser beam moves at a constant velocity, v L along the x-coordinate direction, the time dependent power density can be written as:

Boundary Conditions
All boundaries and surfaces except the top surface where the laser beam impinges are treated as thermally insulated. The boundary conditions for those surfaces and boundaries are given by wheren is the unit normal vector at the boundary. The top surface where laser irradiates is represented with a combination of radiation and convection heat transfer boundary condition as follows: where h c is the convective heat transfer coefficient, T s is the surface temperature, ε is the emissivity of the material, and σ is the Stefan-Boltzmann constant. The choice of the convection coefficient is described in detail in the numerical modeling section.

Numerical Simulation
A three dimensional rectangular slab of PMMA was considered to perform numerical simulation in COMSOL 5.2a as shown in Figure 1a. The dimensions of the rectangular slabs were parametrically defined based on the size of the laser beam and channel length to ensure that the region affected by heat is entirely contained within the model.
The time dependent heat equation (Equation (1)), was solved for the domain shown to find the temperature distribution by the laser. The PARDISO, an efficient robust direct solver with backward differentiation formula (BDF) time stepping method, was employed. The temperature distribution was then used to determine where the material was removed. The material was assumed to be removed once the material reached the temperature at which phase change from solid to vapor is complete. Although the phase change and material removal for a polymer such as PMMA will occur over a range of temperatures, the heating is assumed to happen rapidly enough that it can be approximated as a step change between solid and vapor phases at temperature, T = 700 K for PMMA [32].

Numerical Simulation
A three dimensional rectangular slab of PMMA was considered to perform numerical simulation in COMSOL 5.2a as shown in Figure 1a. The dimensions of the rectangular slabs were parametrically defined based on the size of the laser beam and channel length to ensure that the region affected by heat is entirely contained within the model. The time dependent heat equation (Equation 1), was solved for the domain shown to find the temperature distribution by the laser. The PARDISO, an efficient robust direct solver with backward differentiation formula (BDF) time stepping method, was employed. The temperature distribution was then used to determine where the material was removed. The material was assumed to be removed once the material reached the temperature at which phase change from solid to vapor is complete. Although the phase change and material removal for a polymer such as PMMA will occur over a range of temperatures, the heating is assumed to happen rapidly enough that it can be approximated as a step change between solid and vapor phases at temperature, T = 700 K for PMMA [32]. In this study, a convective heat transfer coefficient was used to model phase change and the heat loss associated with the laser ablation. At the melting phase, heat convection occurs at the solid-liquid boundary as well as between the liquid and the surrounding air. Once the liquid begins to vaporize, it boils out, quickly expands, and escapes with large amounts of heat due to convection. This phenomenon resembles pool boiling heat transfer. In pool boiling, boiling occurs on the solid-liquid interface when the solid surface temperature exceeds the saturation temperature of the liquid. The phase transition in pool boiling results in very high levels of convective heat transfer [34]. Therefore, a very high convective heat transfer coefficients on the order of 10 5 W m −1 K −1 was used in order to model laser micromachining in this study. A range of heat transfer coefficients was explored.

Convergence Study
Second order unstructured tetrahedral mesh elements were used to discretize the domain. Convergence analysis was conducted in order to determine the optimal element size to use. Temperature and depth of cut were used as convergence study parameters. Using constant parameters, which are listed in Tables 1 and 2, the location of the T = 700 K isotherm in the domain was recorded at a given time for element sizes ranging from 32 µm (coarse) to 1 µm (very fine) as shown in the Figure 2a. It was found that the profile of the isotherm varied negligibly for element sizes 4 µm and lower. The depth of cut, reported as the largest depth at which temperature reaches removal temperature, was also found to not vary significantly for element sizes 4 µm or lower, as shown in Figure 2b. Therefore, an element size of 4 µm was used for the remainder of the study. In this study, a convective heat transfer coefficient was used to model phase change and the heat loss associated with the laser ablation. At the melting phase, heat convection occurs at the solid-liquid boundary as well as between the liquid and the surrounding air. Once the liquid begins to vaporize, it boils out, quickly expands, and escapes with large amounts of heat due to convection. This phenomenon resembles pool boiling heat transfer. In pool boiling, boiling occurs on the solid-liquid interface when the solid surface temperature exceeds the saturation temperature of the liquid. The phase transition in pool boiling results in very high levels of convective heat transfer [34]. Therefore, a very high convective heat transfer coefficients on the order of 10 5 Wm −1 K −1 was used in order to model laser micromachining in this study. A range of heat transfer coefficients was explored.

Convergence Study
Second order unstructured tetrahedral mesh elements were used to discretize the domain. Convergence analysis was conducted in order to determine the optimal element size to use. Temperature and depth of cut were used as convergence study parameters. Using constant parameters, which are listed in Table 1 and Table 2, the location of the T =700 K isotherm in the domain was recorded at a given time for element sizes ranging from 32 um (coarse) to 1 um (very fine) as shown in the Figure 2a. It was found that the profile of the isotherm varied negligibly for element sizes 4 m and lower. The depth of cut, reported as the largest depth at which temperature reaches removal temperature, was also found to not vary significantly for element sizes 4 m or lower, as shown in Figure 2b. Therefore, an element size of 4 m was used for the remainder of the study.

Property
Units

Results and Discussions
The three dimensional model described in the previous section was used to investigate the effect of laser system parameters, such as laser power, scanning speed, and beam focus on the depth, width, and shape of the channel. The effect of thermophysical properties, such as specific heat, thermal conductivity, and convection heat transfer coefficient on the channel geometry was also investigated. The system parameters and material properties investigated are provided in Tables 1 and 2 respectively. All of the results presented in the report use the constant property values found in these tables unless stated otherwise. The material properties used to correspond to typical acrylic/PMMA properties from literature [32].

Effect of Laser Power and Unfocused Beam
The laser beam is said to be fully focused when the substrate is located at the focal point of the focusing lens. In focused beam laser ablation, the beam converges to its minimum radius value so that the spot size (i.e., the size of the area irradiated by the beam) is minimized. Thus the spot size can be controlled by adjusting the distance between the substrate and focal point using the following equation.
Here, w 0 is the spot size at the focal point. The value of d corresponds to a given spot size, w is shown in Figure 3. The spot size increases as the beam becomes less focused. This becomes less focused when the substrate is further from the focal point as seen in Figure 3. This means that the spot size of the laser beam on the material can effectively be set by adjusting the distance between the lens and the substrate, without needing to adjust the lens itself.
The results of varying levels of focus, starting with a fully focused beam (spot size w = w 0 = 92.5 µm), are shown in Figure 4 for various laser power settings. The results show that as the laser beam becomes less unfocused, channel width increased while depth decreased. This is due to the energy of the laser beam being distributed over a larger area, decreasing the intensity of the heat added at any given point. The decreased intensity prevents the heat from penetrating deeper into the material. For a given spot size, the depth of cut is not perfectly linear with the power. At lower power settings and larger spot sizes, the intensity of the heat was too low to remove material. Therefore there is a minimum laser power needed for creating a channel on the substrate.
Here, w0 is the spot size at the focal point. The value of d corresponds to a given spot size, w is shown in Figure 3. The spot size increases as the beam becomes less focused. This becomes less focused when the substrate is further from the focal point as seen in Figure 3. This means that the spot size of the laser beam on the material can effectively be set by adjusting the distance between the lens and the substrate, without needing to adjust the lens itself. The results of varying levels of focus, starting with a fully focused beam (spot size w = w0 = 92.5 um), are shown in Figure 4 for various laser power settings. The results show that as the laser beam becomes less unfocused, channel width increased while depth decreased. This is due to the energy of the laser beam being distributed over a larger area, decreasing the intensity of the heat added at any given point. The decreased intensity prevents the heat from penetrating deeper into the material. For a given spot size, the depth of cut is not perfectly linear with the power. At lower power settings and larger spot sizes, the intensity of the heat was too low to remove material. Therefore there is a minimum laser power needed for creating a channel on the substrate.  Figure 4b shows cross-sectional profiles of the channels cut by laser at 40% laser power. The larger the spot size, the shallower the channel becomes. As the heat is spread out over a much larger area, increasingly shallow and wide channels are seen. By varying the spot size, either by changing the vertical distance between the laser and the material, or by varying the optical set up, the desired depth and width of channel can be achieved. It is also shown that the channel closely matches the shape of a Gaussian function at smaller spot sizes (92.5 and 150 µm). For smaller spot sizes, the laser beam is more focused with less loss of the impinged power at the surface. However, as the radius increases to 500 µm, the shape deviates from the Gaussian profile, especially at the edge of the top surface. shows cross-sectional profiles of the channels cut by laser at 40% laser power. The larger the spot size, the shallower the channel becomes. As the heat is spread out over a much larger area, increasingly shallow and wide channels are seen. By varying the spot size, either by changing the vertical distance between the laser and the material, or by varying the optical set up, the desired depth and width of channel can be achieved. It is also shown that the channel closely matches the shape of a Gaussian function at smaller spot sizes (92.5 and 150 um). For smaller spot sizes, the laser beam is more focused with less loss of the impinged power at the surface. However, as the radius increases to 500 um, the shape deviates from the Gaussian profile, especially at the edge of the top surface.

Effect of Scanning Speed
During the laser machining, the speed of the laser is generally held fixed for consistency in results. However the scanning velocity of the laser beam during the ablation process directly determines the exposure time of the material. Thus, the laser speed has an impact on the total heat added, and by extension the depth of cut and other characteristics of the microchannel [27]. The results of varying the laser scan speed and fraction of laser power on the depth of cut are presented in Figure 5.

Effect of Scanning Speed
During the laser machining, the speed of the laser is generally held fixed for consistency in results. However the scanning velocity of the laser beam during the ablation process directly determines the exposure time of the material. Thus, the laser speed has an impact on the total heat added, and by extension the depth of cut and other characteristics of the microchannel [27]. The results of varying the laser scan speed and fraction of laser power on the depth of cut are presented in Figure 5.
During the laser machining, the speed of the laser is generally held fixed for consistency in results. However the scanning velocity of the laser beam during the ablation process directly determines the exposure time of the material. Thus, the laser speed has an impact on the total heat added, and by extension the depth of cut and other characteristics of the microchannel [27]. The results of varying the laser scan speed and fraction of laser power on the depth of cut are presented in Figure 5. As the laser velocity was increased, the channel depth decreased. At higher laser scanning speeds, the total amount of energy imparted on the material is reduced, thus resulting in lower cut As the laser velocity was increased, the channel depth decreased. At higher laser scanning speeds, the total amount of energy imparted on the material is reduced, thus resulting in lower cut depths. At this point, it is important to note that the same depth can be found for various power and speed combinations. For example, a power setting of 10% and laser speed of 100 mm/s gives the same depth as setting the power to 20% with a laser speed of 200 mm/s. Laser scanning speed also contributes to the quality of the surface of the channel. For instance, the roughness of the cut has been reported to increase with decreasing laser scanning speeds [27]. The overall cross-sectional profile of the channel also varies with the speed as seen in Figure 5b. It is shown that the width of cut remains approximately the same, regardless of the scanning speed, while the depth varies. Thus, various aspect ratios of channel cross-section can be obtained simply by varying the laser speed.

Effect of Thermal Conductivity
The effect of substrate thermal conductivity on the depth of cut is presented in Figure 6. Thermal conductivity values for polymers usually fall in the range of 0.1-0.3 W m −1 K −1 [37]. For instance, the thermal conductivity of PMMA is 0.19 W m −1 K −1 . In laser ablation, heat is localized to a very small region of intense heat. In a material with poor thermal conductivity, heat is not able to spread far from the initially irradiated region. Figure 6 shows that higher thermal conductivity results in a deeper cut. Heat can easily spread further down into the material with higher thermal conductivity before it is ultimately removed by convection. depths. At this point, it is important to note that the same depth can be found for various power and speed combinations. For example, a power setting of 10% and laser speed of 100 mm/s gives the same depth as setting the power to 20% with a laser speed of 200 mm/s. Laser scanning speed also contributes to the quality of the surface of the channel. For instance, the roughness of the cut has been reported to increase with decreasing laser scanning speeds [27]. The overall cross-sectional profile of the channel also varies with the speed as seen in Figure 5b. It is shown that the width of cut remains approximately the same, regardless of the scanning speed, while the depth varies. Thus, various aspect ratios of channel cross-section can be obtained simply by varying the laser speed.

Effect of Thermal Conductivity
The effect of substrate thermal conductivity on the depth of cut is presented in Figure 6. Thermal conductivity values for polymers usually fall in the range of 0.1-0.3 Wm −1 K −1 [37]. For instance, the thermal conductivity of PMMA is 0.19 Wm −1 K −1 . In laser ablation, heat is localized to a very small region of intense heat. In a material with poor thermal conductivity, heat is not able to spread far from the initially irradiated region. Figure 6 shows that higher thermal conductivity results in a deeper cut. Heat can easily spread further down into the material with higher thermal conductivity before it is ultimately removed by convection.

Effect of Specific Heat
The effect of the specific heat of substrates on the depth of cut was investigated for a range of constant specific heats (between 500-3000 Jkg −1 K −1 ) as well as temperature dependent specific heat. Specific heat varies significantly for a material that undergoes a phase change. Temperature dependent specific heat was specified using the following equation adopted from [32].

Effect of Specific Heat
The effect of the specific heat of substrates on the depth of cut was investigated for a range of constant specific heats (between 500-3000 J kg −1 K −1 ) as well as temperature dependent specific heat. Specific heat varies significantly for a material that undergoes a phase change. Temperature dependent specific heat was specified using the following equation adopted from [32].
where H(T) is the smoothed Heaviside function centered at the vaporization temperature, and D(T) is a function used to distribute the latent heat of phase change across the phase change interval. Smoothed step functions were used when phase change occurs. Figure 7 shows that for constant lower specific heat values, the depth is much more sensitive to the effect of laser power. At higher laser powers, the effect is more pronounced. It can be noted that for a range of depths, multiple combinations of specific heat and laser powers are available. For instance, a material of specific heat 500 J/kg·K with a laser power of 6.5 W gives the same channel depth as a using 19.5 W laser power on a material with specific heat 1500 J/kg·K. This means that the same depth of cut is achievable across various types of substrate materials through adjusting the laser power. The depth of cut resulting from the phase dependent specific heat also shown in Figure 7. It shows that the constant specific heat assumption can overestimate the depth of cut by as much as 10%. This is significant for microfluidic applications, so it is imperative to consider temperature dependent specific heat in computer simulation of laser micromachining of microfluidic devices.
Micromachines 2019, 10, x 9 of 12 same depth of cut is achievable across various types of substrate materials through adjusting the laser power. The depth of cut resulting from the phase dependent specific heat also shown in Figure 7. It shows that the constant specific heat assumption can overestimate the depth of cut by as much as 10%. This is significant for microfluidic applications, so it is imperative to consider temperature dependent specific heat in computer simulation of laser micromachining of microfluidic devices.

Effect of Convective Heat Transfer Coefficients
The effect of convective heat transfer coefficients on the depth of cut was evaluated for a wider range of convection coefficients between 250 and 5 × 10 5 Wm −2 K −1 . Typical forced convection heat transfer coefficient of 250 Wm −2 K −1 was used to compare with the significantly higher values. Figure 8 shows that use of 250 Wm −2 K −1 as heat convection coefficients cannot capture the variation of the channel for various laser power. The variation of depth due to increasing laser power becomes even more insignificant when lower convective coefficients were used (not shown). However, it has already been established previously that laser power has a significant effect on the depth of cut in laser ablation [21,25,29]. The depth of cut decreases with increased convection coefficient, regardless of the power settings. However, the depth of cut increases significantly with laser power at higher convection coefficients as well, which resembles the experimental observation of laser ablation. The larger heat convection coefficient accounts for not only the amount of heat removed by natural convection at the surface, but also the effect of additional heat escaping the bulk material. Heat escapes during melting and vaporization of the material due to laser treatment.

Effect of Convective Heat Transfer Coefficients
The effect of convective heat transfer coefficients on the depth of cut was evaluated for a wider range of convection coefficients between 250 and 5 × 10 5 W m −2 K −1 . Typical forced convection heat transfer coefficient of 250 W m −2 K −1 was used to compare with the significantly higher values. Figure 8 shows that use of 250 W m −2 K −1 as heat convection coefficients cannot capture the variation of the channel for various laser power. The variation of depth due to increasing laser power becomes even more insignificant when lower convective coefficients were used (not shown). However, it has already been established previously that laser power has a significant effect on the depth of cut in laser ablation [21,25,29]. The depth of cut decreases with increased convection coefficient, regardless of the power settings. However, the depth of cut increases significantly with laser power at higher convection coefficients as well, which resembles the experimental observation of laser ablation. The larger heat convection coefficient accounts for not only the amount of heat removed by natural convection at the surface, but also the effect of additional heat escaping the bulk material. Heat escapes during melting and vaporization of the material due to laser treatment.

Conclusions
This paper presents a mathematical model where high convective heat transfer coefficients were used to simplify the modeling of complex phase transitions and material removal processes found in laser microfabrication. The developed model was used to determine the effect of individual laser system parameters and material properties on the resulting fabricated channel. It was found that the laser system parameters such as power, beam focus, and scanning speed can be used to control and, eventually fine-tune the fabrication process for both channel width and height. The laser impingement spot size can be controlled by the level of focus of the beam. Similarly by selecting appropriate substrate materials, it would be possible to fabricate a channel with a desired width and depth using laser machining. This study provides the following conclusions: 1. The laser impingement spot size depends on the level of focus of the beam. Smallest spot sizes are achieved with a fully focused beam when material is located at its focal point, and spot size can be increased by increasing the distance between laser and material.
2. Larger laser spot size creates wider and shallower channels due to more scattered nature of the impinged energy.
3. The laser scanning speed mainly affects the depth with little to no impact on the width of the channel. If the laser scanning speed is increased, the depth of cut is decreased.
4. Substrates with higher thermal conductivity will have a deeper depth of cut for a given laser power and other system parameters.
5. The depth of cut increases as the laser power increases. Larger heat convection coefficient can be used to simplify the mathematical modeling of laser machining, to predict qualitatively the variation in depth as a function of laser power. 6. The higher the specific heat of a substrate the lower the depth of cut for a given laser power. The use of temperature dependent specific heat is more appropriate for computer modeling compared to constant specific heat in laser micromachining.