Selective laser foaming is a novel process that combines solid-state foaming and laser ablation to fabricate an array of microliter tissue engineering scaffolds on a polymeric chip for biomedical applications. In this study, a finite element analysis (FEA) model is developed to investigate the effect of laser processing parameters. Experimental results with biodegradable polylactic acid (PLA) were used for validation. It is found that foaming always occurs before ablation, and once it occurs, the temperature increases dramatically due to an enhanced laser absorption effect of the porous structure. The geometry of the fabricated scaffolds can be controlled by laser parameters. While the depth of scaffolds can be controlled by laser power and lasing time, the diameter is more effectively controlled by the laser power. The model developed in this study can be used to optimize and control the selective foaming process.

## Introduction

Selective laser foaming is a novel material processing technique developed to fabricate an array of microliter tissue engineering scaffolds for high-throughput drug screening tests [1–3]. Two-dimensional (2D) cell cultures are widely used in the drug discovery and development process. However, cells growing on 2D substrates are forced to form an artificial monolayer, losing their in vivo three-dimensional (3D) morphology and the ability to produce realistic drug testing results [4,5]. It is well known that for cells to recapitulate tissue specific functions, a tissue engineering scaffold is needed to support the reformulation of 3D tissue structure [6,7]. Given the fact that experimental drugs are often available in minute amounts, these tissue engineering scaffolds need to be small and fabricated in a large array such that high-throughput tests could be performed in parallel to reduce the time and cost, as well as the statistical variability, of the drug development process.

Figure 1 shows a schematic of the developed selective laser foaming process [1]. It combines the characteristics of laser ablation and solid state foaming, where a polymer sample is first saturated in a high pressure vessel containing an inert gas such as nitrogen (N_{2}) and carbon dioxide (CO_{2}). Under a certain saturation condition, the sample is retrieved from the pressure vessel and treated with laser irradiation. The heat induced by laser softens the polymer matrix and lowers the solubility of the inert gas, causing gas bubbles to nucleate and grow and the formation of porous structure. As the heat is continuously generated, the center of foamed region undergoes photothermal ablation, creating an inverse cone-shaped porous well that is conducive for cell culturing. The selective laser foaming process has been demonstrated with both poly(methyl methacrylate) and poly(lactic acid) (PLA) [1,2]. The volume of scaffolds fabricated was on the level of a few microliters. Human glioblastoma cells (T98) were successfully cultured in these scaffolds, forming organotypic clusters that could not be obtained in 2D cell culture.

In this study, a finite element analysis (FEA) model is developed to evaluate parameter effects of the selective laser foaming process. The ultimate goal is to control the size and consistency of fabricated tissue engineering scaffolds, such that the process could be used with a range of materials and various cell types that may require different geometric features of the scaffolds. With the FEA model, the effects of laser power and lasing time on the size and shape of the fabricated scaffolds were investigated.

## Finite Element Modeling

Finite element modeling is frequently performed to study complex phenomena in laser-based manufacturing processes [8–10]. However, little has been done on modeling of an integrated laser foaming and ablation process. Figure 2 shows the model geometry and boundary conditions of the FEA model developed in this study. The polymer sample and laser beam are axisymmetric along the center of the beam. Therefore, the model is simplified into a 2D model, and only half of the domain is considered.

where *h* is a convective heat transfer coefficient, *T*_{surf} is the temperature at the sample surface, *T*_{0} is the ambient temperature, and $n$ is the normal vector of the corresponding boundary. The adiabatic condition is imposed on boundary 4. Due to the effect of the laser cooling fan, a typically forced convective heat transfer coefficient for moderate speed air flow, 100 W/(m^{2} K) [11], was used in the simulation.

*r*and

*z*are the spatial coordinates and

*ρ*,

*C*, and

_{p}*k*are density, specific heat, and thermal conductivity of the material, respectively. The heat source

*Q*(

*r, z*) is the energy input induced by the laser. It is assumed that the laser intensity in the radial direction follows a Gaussian distribution and the absorption of laser energy by the polymer can be defined by the Beer–Lambert law [13]

where *α* and *β* are the absorption coefficient and the attenuation coefficient of the polymer, and *w* is the laser beam width.

_{00}beam, since the target substrate is placed within the Rayleigh range. The divergence is modeled using the following equation [14]:

where *x* is the distance between the laser outlet and the sample, *w*(*x*) is the beam diameter at the distance of *x*, *w*_{0} is the beam diameter at the laser outlet, and *λ* is the laser wavelength.

*ρ*) of the polymer is reduced. At the same time, the absorption coefficient (

*α*) and attenuation coefficient (

*β*) are increased, because the bubbles scatter the laser and more laser energy is absorbed by the polymer. The material properties after the phase changes including the density (

*ρ*), absorption coefficient (

*α*), and attenuation coefficient (

*β*) are defined as follows:

where *χ* denotes any of the above three material properties, *χ _{s}* is the material property of unfoamed polymer,

*χ*is that of foamed polymer,

_{f}*χ*is that after ablation,

_{A}*T*and

_{f}*T*are the starting temperatures of foaming and ablation [15,16], respectively, and

_{A}*ΔT*

_{1}and

*ΔT*

_{2}are the temperature ranges of foaming and ablation, respectively. Smoothed Heaviside step functions between phases are used to avoid extreme discontinuity and reduce convergence error in simulation.

*H*is the latent heat of phase change,

*T*

_{start}and

*T*

_{end}are the starting and ending temperatures of the phase change, respectively. Radice et al. [17] and Groulx and Ogoh [18] suggested to define one parameter that can be used to track the progress of material over phase change. The specific heat of the polymer is thus defined over different temperature ranges as

where *C*_{ps} is the specific heat of unfoamed polymer, *C*_{pf} is the specific heat of foamed polymer, *C*_{pA} is the specific heat of ablated polymer, which has a similar value to air, *C _{p}*

_{1}is the changing specific heat during foaming, which follows a smoothed Heaviside step function from

*C*

_{ps}to

*C*

_{pf},

*C*

_{p}_{2}is the changing specific heat during ablation that can be expressed again as a smoothed Heaviside step function from

*C*

_{pf}to

*C*

_{pA}, and

*d/dT*represents a Dirac delta function which is a derivative of the smoothed Heaviside step function.

In this study, CO_{2} was used as the blowing agent and biodegradable PLA as the polymer matrix material. The latent heat of foaming was calculated using the latent heat of CO_{2} sublimation and gas concentration of the saturated polymer sample. The latent heat of CO_{2} sublimation was chosen to be 25.2 kJ/mol [19]. The gas concentration was determined based on previous studies of solid state foaming of PLA [20–22]. To account for bubble growth in the PLA matrix, it was assumed that a certain amount of sublimation energy was consumed by CO_{2} to expand the polymer matrix. This amount was determined to be 10% by matching the modeling and experimental results. The latent heat of ablation was obtained from a published thermal degradation study of PLA [23]. The parameters used in the model are summarized in Table 1.

The FEA model was implemented using the commercial software COMSOL Multiphysics^{®} (COMSOL, Burlington, MA) with a backward differentiation and automatic time stepping procedure for computation. An extremely fine mesh was used, and the minimum element size was 50 *μ*m. The heat source was applied using a user-defined function of location and temperature. Once the entire computation process was completed, the temperature distribution was stored and utilized to estimate the foamed and laser-ablated regions within the simulation domain. The temperature of each element was compared to the foaming and ablation temperatures of PLA to determine the foaming and ablation regions. Figure 3 shows a block diagram of the solution procedure.

## Experiments

The model developed in this study was validated with experimental data. Substrates used in the experiment were disk samples having a diameter of 30 mm and average thickness of 1.5 mm. The samples were injection-molded with PLA powder (ECORENE^{™} NW40) which is a mixture of D/L product from ICO Polymers (Allentown, PA). The PLA samples were placed in a high pressure CO_{2} chamber until they reach the fully saturation condition. The saturation pressure was 2 MPa, and the averaged concentration of CO_{2} in the PLA samples was 8.5 wt % after the saturation process. After they were retrieved from the pressure vessel, the gas-impregnated samples were exposed under continuous laser irradiation by a CO_{2} laser (Synrad Firestar v30 from Synrad, Inc., Mukilteo, WA). For comparison, unsaturated samples were treated with the same laser conditions, with which only ablation was generated. Table 2 summarizes the sample conditions and laser parameters used in the experiment. The treated samples were freeze-fractured in liquid nitrogen. A scanning electron microscope (SEM) was used to obtain cross-sectional images of the treated samples.

## Results and Discussion

### Model Validation.

Laser attenuation in PLA is not well defined in literature due to its various crystallinity and L/D composition [24]. Therefore, the attenuation coefficients of both unfoamed and foamed PLA were estimated with the ablation and foaming results with a laser power of 7.7 W and lasing time of 1 s. The attenuation coefficients were varied in the simulation, and the predicted ablation and foaming regions were compared with the cross-sectional SEM images of the samples. The ablation and foaming regions were identified based on the ablation and foaming temperatures of PLA. Until the results were matched, these attenuation coefficients were determined to be 200 mm^{−1} for solid PLA and 1075 mm^{−1} for foamed PLA. Figure 4 shows the comparison between experimental and simulation results for both unfoamed and foamed samples under these conditions. The laser-ablation profiles were determined by the thermal degradation temperature of PLA. Figure 4(d) shows the temperature distribution of foamed region predicted using the solid state foaming temperature of PLA. In Figs. 4(c) and 4(d), the elements in the ablated region were removed from the model.

Once determined, the attenuation coefficients were used for the rest laser foaming conditions, with the laser power ranging from 2.0 W to 10.3 W and lasing time from 0.25 s to 1.00 s. Figure 5 shows the comparison between the modeling and experimental results. The predicted ablation and foaming profiles were overlaid with the corresponding cross-sectional SEM images from experiments. For all the laser foaming conditions, the predicted ablation and foaming profiles showed good agreement with the experimental results.

In order to examine the combined effects of laser power and lasing time, the relationship between the volume of foamed region and laser energy is plotted in Fig. 6 for both the simulation and experimental data. The volume of the foamed region was considered to be cone shaped. The prediction from the FEA model is in good agreement with the experimental data at low laser energy. As the laser energy increases, the volume of the foamed region is affected by the sample thickness, and therefore may not be accurately predicted. This however will not affect the applicability of the FEA model, since high laser energy settings will generate a through hole in the sample, a situation which should be avoided in normal operations of the process.

### Predicted Temperature Distribution.

Predicted temperature distribution at the sample surface is shown in Fig. 7. The laser power used in the simulation was 7.7 W. It is seen that the temperature profile changes significantly with time. Since the laser power density has a Gaussian distribution, the temperature at the center is much higher than that in the outer region. In addition, two phase change regions are observed for foaming and laser ablation, respectively. The first phase change happens at about 340 K when foaming starts. The second happens around 440 K when PLA starts to decompose. During the two phase changes, the temperature at the center of the laser irradiated region continues to rise. After foaming, the center temperature increases dramatically due to enhanced laser absorption by the foamed structure.

### Effect of Laser Power and Lasing Time on the Size of Ablated Region.

The simulation model was used to study the effects of laser power and lasing time on the ablated region of foamed samples. The depth of ablated region as a function of lasing time is shown in Fig. 8. Before it reaches the substrate thickness limit, the ablation depth has a linear relationship with lasing time. The slope of this linear trend is determined by the laser power. If the laser power is higher, the depth of ablated region reaches its limit faster. For example, a through-hole can be created within 0.8 s at a 10.3 W laser power (16% of the max laser power), while a 7.7 W laser power (12% of the max laser power) will not be able to create a through-hole even after 1 s of lasing time. For a lower power, e.g., 4.9 W (8% of max laser power), ablation will not start until 0.1 s after the onset irradiation. For a 2.0 W laser power (4% of max laser power), there is no ablation. The diameter of the ablated region as a function of lasing time is shown in Fig. 9. Similar trends have been observed. The diameters of ablated regions approach to a certain limit. However, this limit is not caused by the thickness, but different power settings of the laser.

### Effect of Laser Power and Lasing Time on the Size of Foamed Region.

The simulation model was also used to study the effect of laser power and lasing time on the size of foamed regions. The depths of the foamed region are plotted with lasing time in Fig. 10, which shows similar characteristics as the ablation depth seen in Fig. 8. However, there is almost no delay in foaming. Foaming occurs as soon as the laser starts to irradiate on the samples. This is because the temperature threshold for foaming is much lower than that for ablation.

The diameters of foamed regions under different power settings are shown in Fig. 11. The foamed diameters are plotted in terms of lasing time, and each line represents a different laser power. Similar to the ablated diameter, the relationship between the foamed region diameter and lasing time follows a logistic function. However, the diameter of foamed region is much larger and approaches to its steady state value much faster. For example, the 7.7 W laser power generates a foamed region about 4 mm in diameter, whereas the ablation region diameter is only 2.5 mm. In addition, the final diameters of foamed regions under different laser powers are also different, similar to the trend seen in Fig. 9. This is due to the Gaussian distribution of laser power density, which causes a different spot size at any given laser power setting.

As seen in Figs. 9 and 11, both foaming and ablation occur rapidly in the initial stage of the process, and the diameter approaches to a certain equilibrium value at a given laser power. A higher laser power tends to yield a larger foamed and ablated region. This is because that the diameters of the foamed and ablated regions are affected by the distribution of the laser power density. As shown in Fig. 12, assuming the same power density will lead to the increase in same temperature, the foaming and ablation temperatures will impose two different thresholds on the Gaussian distribution of the laser power density. Different power settings will lead to different diameters of the foamed and ablated regions. Therefore, laser power instead of lasing time should be used to obtain scaffolds with different diameters.

## Conclusions

A finite element model has been developed to study the selective laser foaming process. The model was validated with experimental data. It has been found that once foaming occurs, the temperature increases dramatically because of the enhanced laser absorption effect of the porous structure. For a fixed laser power, the depths of both ablated and foamed regions are linearly proportional to the lasing time before they reach a limit defined by the substrate thickness. Laser power has an effect on the slope of this linear trend. The diameters of the ablated and foamed regions increase rapidly at the early stage of the process and approach to a limit determined by the laser power setting. Therefore, laser power is a more effective variable to use in order to control the diameter of the foamed and ablated regions. Laser foaming is a complex process involving heating, foaming, and ablation. The FEA model developed in this study provides a useful tool to aid future development of the process.

## Funding Data

National Science Foundation (Grant No. CMMI-1131710).