Open Access
3 May 2017 Direct numerical simulation of the initial stage of a thermally induced microcavitation in a water-rich biotissue triggered by a nanosecond pulsed laser
Sy-Bor Wen, Kevin J. Ly, Arun Bhaskar, Morgan S. Schmidt, Robert J. Thomas
Author Affiliations +
Abstract
A numerical analysis capable of describing the early stage of a thermal microcavitation process in a water-rich biotissue without avalanche breakdown was developed. The analysis successfully reproduced the laser-induced heating, vapor bubble formation, bubble expansion, and shockwave propagation inside a water-rich biotissue during a thermal microcavitation process. Based on the analysis, it was determined that the evolution of the temperature, pressure, and laser-induced shockwave is dependent on the incident laser energy and laser pulse width. On the other hand, the early stage dynamics of the microcavitation process showed little dependence on the elastic modulus of the biotissue for the laser and tissue conditions studied.

1.

Introduction

Laser-induced microcavitation refers to the rapid formation and expansion of a vapor bubble inside the biotissue when it is exposed to intense enough pulsed laser energy. The vapor bubble formation can be purely thermal (i.e., due to light absorption of the tissue in condensed or vaporized phases) with or without avalanche breakdown. With the associated microscale dissection left in the biotissue, laser-induced microcavitation is a common approach for high-precision biosurgeries. For example, laser-induced microcavitation is used in laser in situ keratomileusis (LASIK) to precisely reshape the midstromal tissue of the cornea using an excimer laser beam. Several studies were conducted since the first LASIK surgery in order to understand the laser-induced microcavitation process. For example, Vogel et al.13 studied the formation, oscillation, and collapse of vapor bubbles at time steps after a nanosecond (ns) and femtosecond (fs) laser pulse during a microcavitation process. Tinne et al.4 studied the interaction of the spatially separated vapor bubbles from a microcavitation event. Blake and Gibson5 investigated the dynamics of cavitation bubbles near boundaries. Kodama and Tomita6 studied bubble–shockwave interaction during a microcavitation process. Schmidt et al. investigated the microcavitation processes of retinal pigment epithelium melanosomes and polystyrene microspheres.79 Most of these studies were based upon experimental methods and relied on fast laser imaging to capture the transient dynamics of the vapor bubble(s) under different laser, tissue, and event boundary conditions.

Among the few theoretical analyses that explore the laser-induced microcavitation process, we find most focus on either the laser-induced heating of the biotissue before the microcavitation threshold or the oscillation of the vapor bubble inside the biotissue that occurs a relatively long time (>10  s of μs) after the occurrence of microcavitation.1013 Creating models for microcavitation processes would be beneficial in areas of heating and stress during microcavitation.14,15 The subsequent modeling and identification of these physical events would help optimize use conditions for laser surgery associated with damage promoting events such as reduction of heat zones and lessening the peak stress values during laser-induced microcavitation.16 In this study, a direct simulation model of an early microcavitation process due to laser absorption with axial geometric symmetry is performed in water-rich biotissue. Using the early stage analysis, we were able to determine the heat and stress affected zones during a microcavitation, which is important in the optimization of laser surgery. Furthermore, with the capability to relate the temperature, pressure, phase, and stress fields inside the biotissue as functions of physical properties of the tissue and the laser conditions, the theoretical analysis can serve as a tool to identify the physical properties of the biotissue from the measured temperature, pressure, phase, and stress fields during a microcavitation. This reverse approach can be valuable, especially when physical properties of certain biostructures (e.g., the absorption coefficient of cell organelles) are difficult to measure directly with other methods.

In this study, we focus on the thermal microcavitation process without plasma formation and avalanche breakdown. In other words, the microcavitation event occurs by strong optical absorption in the biotissue and the followed rapid phase change with a maximum temperature less than 1500 to 2000 K. This specific type of thermal microcavitation by long wavelength nanosecond lasers (e.g., Er-YAG laser) occurs in water-rich biotissues, commonly applied in laser surgeries. Compared with other microcavitation processes, this specific type of thermal microcavitation requires less laser energy and shows high consistency with smoother dissection structure during the experiment.17 To model the occurrence of this specific thermal microcavitation in biotissues, mass, momentum, and energy conservation, equations are required.1,15 In addition, pressure (P) and specific heat (cv) as functions of temperature (T) and density (ρ, which is inverse of specific volume) of the biotissue should be specified, which can be challenging for a large temperature variation simulation since both of them can vary significantly and nonlinearly with respect to temperature and density. To solve this issue, we adopted a cubic form equation of state of water as a first-order approximation of the P, ρ, T correlation of the biotissue when it is in liquid and vapor phases. When the tissue is in an elastic solid phase (for lower temperature), a scaling factor is used to adjust the equation of state to ensure the analytical expression provides a bulk modulus (i.e., inverse of compressibility) similar to the experimental value. Once the P, ρ, T correlation is specified, the specific heat cv of the water-rich biotissue can then be estimated with Maxwell’s thermodynamic relations (Sec. 2.2).18 The description of the new theoretical analysis capable of describing the occurrence of microcavitation in water-rich biotissue during a thermal microcavitation process without plasma formation and avalanche breakdown is in Sec. 2. Section 3 includes the obtained results and discussion.

2.

Theoretical Analysis

In this section, we first present the proposed equation of state for water-rich biotissues as well as the associated analytical expression of the internal energy as a function of temperature and specific volume in each phase. We then present the required governing equations and constitutive equations to solve the temperature, pressure, specific volume, and stress fields inside the biotissue during the occurrence of microcavitation.

2.1.

Equation of State for Water-Rich Biotissues

In this study, as a first-order approximation, we assume the equation of the state of water-rich biotissues (e.g., 78% of water in corneal stroma) can be approximated with that of water when the biotissue is in the liquid or vapor phases.15 When the biotissue is in the elastic solid phase under a lower temperature (i.e., less than dissociation temperature of the tissue) and higher density (i.e., around its intrinsic density), the equation of the state of water provides orders of magnitude higher prediction of the pressure value. Therefore, in order to relate the tissue pressure to its temperature and density following the experimentally measured elastic modulus, a scaling factor is required for this particular phase. To illustrate this idea, we start from the modified Peng–Robinson equation of state (PR EOS), which provides a good estimation of the pressure, temperature, and specific volume of pure water at the subcooled liquid and superheated vapor phase.19 This modified PR EOS is adapted to predict the local pressure of water-rich biotissues in this study when the tissue is in liquid and vapor phases. An advantage of using analytical EOS is the continuous derivative values at each state, which can significantly enhance the numerical stabilities. A full expression of the PR EOS is

Eq. (1)

P=RTv˜b˜a(T)v˜2+2b˜v˜b˜2,a(T)=0.45724(RTc)2Pc[1+m(1Tr)]2;m=0.37464+1.5422ϖ0.26992ϖ2,Tr=TTc;b˜=0.0778RTcPc;v˜=v+c,
where P represents the thermal pressure, T is the temperature, v is the specific volume, Tc is the critical temperature of water (647.1  K), Pc is the critical pressure of water (22,060  kPa), ϖ is the acentric factor of water (0.334), c is the corrected factor for specific volume (0.24091×103  m3/kg), and R is the gas constant of water (0.4615  kJ/kg·K). Through Maxwell’s equal area rule, the specific volumes of saturated liquid and vapor as well as the saturated pressure at each temperature can be determined for PR-EOS as shown in Fig. 1. For specific volume within the saturate curve, phase change occurs and the vapor qualify is estimated from the specific volume of the local state with respect to the specific volumes of the saturated liquid and vapor.

Fig. 1

(a) Saturated curve obtained from the modified PE EOS. Left light gray line shows the specific volume of saturated liquid vf. Right dark gray line shows the specific volume of saturated vapor vg. (b) Saturated pressure Psat obtained from the modified PE EOS with Maxwell equal area rule.

JBO_22_5_056002_f001.png

When the biotissue is in the saturated phase (i.e., the specific volume of the biotissue is within vf and vg as shown in Fig. 1), the pressure of the biotissue becomes Psat(T) determined from the modified PR EOS as shown in Fig. 1(b). Note that metastable state during the phase change is not considered in this study considering the ns (and longer time) duration focused in this study allows material to relax to local thermodynamics equilibrium.

When the temperature of the water-rich biotissue is less than its dissociation temperature (e.g., 340  K for stroma of cornea) and the specific volume of the biotissue is still around its intrinsic value (0.001  m3/kg), the fibroblast cells (or collagen in each cell) of the biotissue hold the structure of the tissue as an elastic solid. The bulk modulus B of biotissue in elastic solid phase varies from MPa level (e.g., stroma of cornea) to GPa level (e.g., brainstem tissue), which can be very different from that determined with the modified PR EOS for pure water (i.e., 2  GPa). In other words, the pressure determined with Eq. (1) can have several orders of magnitude of overestimation when the biotissue is in elastic solid phase. To rectify this error, when T is less than the dissociation temperature and when the specific volume of the biomaterial is less than vf as shown in Fig. 1(a), the following equation of state is adopted

Eq. (2)

Pelastic=[RTv˜b˜a(T)v˜2+2b˜v˜b˜2Psat(T)]*(scaling factor)+Psat(T),
where the scaling factor scales down the pressure increment when the specific volume of the tissue is less vf. The scaling factor is determined from the experimentally measured bulk modulus B of the biotissue, which is a function of elastic modulus E and Poisson ratio ν, as

Eq. (3)

(scaling factor)×[RTatmv˜atmb˜a(Tatm)v˜atm2+2b˜v˜atmb˜2Psat(Tatm)]=Bρatmρl(Tatm)ρl(Tatm)=E1+νν12νρatmρl(Tatm)ρl(Tatm),
with Tatm is the room temperature; v˜atm=vatm+c and vatm is the specific volume of the tissue under room condition; Psat(Tatm) is the saturated water pressure at room temperature; ρatm is the density of the tissue under room temperature; ρl(Tatm) is the density of saturated liquid water under room temperature. The scaling factor becomes one when the temperature of the biotissue is higher than the dissociation temperature and when the specific volume of the tissue is still less than vf.

2.2.

Internal Energy (e)

Internal energy e of water-rich biotissue is another property that should be specified as a function of temperature and specific volume before a large temperature variation simulation. In a traditional thermo-optical simulation of biotissue, e is a function of T only and the time variation of e can be expressed as

Eq. (4)

decv(T)dT,
where cv represents the specific heat of the biotissue. However, such correlation is no longer valid when a significant change of density (i.e., inverse of specific volume) of the tissue occurs during a thermos-optical simulation. For the more general situation during a microcavitation process, consider dependence of e with the specific volume with

Eq. (5)

de=cv(T)dT+(ev)Tdv.
With the help of Maxwell identity, (e/v)T can be expressed as a function of temperature and pressure

Eq. (6)

(ev)T=T(PT)VP.
When the biotissue is in the vapor phase or when the temperature of the biotissue is higher than the critical temperature Tc, modified PR EOS can be substituted into Eq. (5) to have an explicit expression of (e/v)T as

Eq. (7)

(ev)T,vapor=Ta(T)/Tv˜2+2b˜v˜b˜2+a(T)v˜2+2b˜v˜b˜2=0.45724(RTc)2Pc[1+m(1Tr)](1+m)v˜2+2b˜v˜b˜2.
Integrating Eq. (7) with respect to v provides the following explicit expression of e valid for v>vg

Eq. (8)

e(T,v)vapor=0.45724(RTc)2Pc[1+m(1Tr)](1+m){ln[b+v˜2bb+v˜+2b]22bln[b+vR+c2bb+vR+c+2b]22b}+fct(T).
With vR selected as a very large specific volume (here, we selected 200  m3/kg) for which the internal energy can be approximated with a pure function of temperature fct(T), the values of fct(T) are obtained from IAPWS-95 equation.20 The associated de in the supercritical region has

Eq. (9)

de(T,v)vapor={0.45724(RTc)2Pc[12m1Tr1Tc](1+m){ln[b+v˜2bb+v˜+2b]22bln[b+vR+c2bb+vR+c+2b]22b}+cR(T)}dT+0.45724(RTc)2Pc[1+m(1Tr)](1+m)v˜2+2b˜v˜b˜2dv,
where cR(T) is the specific heat of the reference high specific volume state, which can also be extracted from the IAPWS-95 equation. When the biotissue is in the saturated phase (i.e., v is between vf and vg of Fig. 1), P=Psat(T) and Eq. (6) becomes

Eq. (10)

(ev)T=[TdPsat(T)dTPsat(T)].
Therefore, by integrating (ev)T from v to vg, we have

Eq. (11)

e(T,v)saturated=[TdPsat(T)dTPsat(T)]·[vvg(T)]+eg(T),
where eg(T) is the internal energy of the saturated vapor at temperature T and specific volume vg. The associated de in the saturated region has

Eq. (12)

de(T,v)saturated={Td2Psat(T)dT2[vvg(T)][TdPsat(T)dTPsat(T)]dvg(T)dT}dT+[TdPsat(T)dTPsat(T)]dv+deg(T),
deg(T) of Eq. (12) can be determined from Eq. (8) with v=vg as

Eq. (13)

deg(T)={{0.45724(RTc)2Pc[12m1Tr1Tc](1+m){ln[b+v˜g2bb+v˜g+2b]22bln[b+vR+c2bb+vR+c+2b]22b}+cR(T)}+0.45724(RTc)2Pc[1+m(1Tr)](1+m)v˜g2+2b˜v˜gb˜2dvg(T)dT}dT.

When the biotissue is in the liquid phase, by integrating Eq. (6) from v to vf, we have

Eq. (14)

e(T,v)liquid=0.45724(RTc)2Pc[1+m(1Tr)](1+m){ln[b+v˜2bb+v˜+2b]22bln[b+vf+c2bb+vf+c+2b]22b}+ef(T),
where ef(T) is the internal energy of the saturated liquid at temperature T and specific volume vf. The associated de in the liquid region has

Eq. (15)

de(T,v)liquid={0.45724(RTc)2Pc[12m1Tr1Tc](1+m){ln[b+v˜2bb+v˜+2b]22bln[b+vf+c2bb+vf+c+2b]22b}0.45724(RTc)2Pc[1+m(1Tr)](1+m)v˜f2+2b˜v˜fb˜2dvfdT}dT+0.45724(RTc)2Pc[1+m(1Tr)](1+m)v˜2+2b˜v˜b˜2dv+def(T),
where def(T) of Eq. (15) can be determined from Eq. (12) with v=vf.

When the temperature is less than the dissociation temperature and the specific volume is higher than vf, water-rich biotissue is in an elastic solid phase and is almost incompressible (i.e., Poisson ratio close to 0.5). As a result, Eq. (5) is reduced to Eq. (4) and the internal energy is almost independent of the specific volume. Therefore, de(T,v)elastic solidde(T,vf)saturated with de(T,vf)saturated from Eq. (12) under the condition v=vf. Note that de(T,v)liquid predicted in Eq. (15) is also similar to de(T,vf)saturated from Eq. (12) when the liquid temperature is low and the compressibility of liquid is close to zero based on the modified PE EOS. In other words, we can use either de(T,v)liquid or de(T,vf)saturated to estimate the values of de(T,v)elastic solid at a low temperature with minimum error.

2.3.

Mass Conversation Equation for the Water-Rich Biotissue

The density-based method adopted in this study better handles the large density variation of the biotissue due to the presence of the vapor bubble during a microcavitation. In order to determine the density evolution of the biotissue, the mass conservation equation [Eq. (16)], with ρ is the local density and V is the local velocity field of the biotissue, is used

Eq. (16)

ρt=·(ρV).

2.4.

Momentum Conservation Equation for the Water-Rich Biotissue

The Navier–Stokes equations and the definition of material derivative in Euler [found in Eq. (17) and (18)] describe the evolution of the velocity field V and the associated displacement field u in the biotissue during the microcavitation process

Eq. (17)

ρVt+ρ(V·)V=·σ+Fb,

Eq. (18)

ut+(V·)u=V.
The gravitation acceleration Fb is neglected in the simulation considering it is not important during the first few microseconds of the highly transient microcavitation process. The stress tensor σ experienced by the biotissue during the microcavitation depends on the phase of the biotissue. When the biotissue is in the elastic material phase, the stress tensor σik can be estimated with linear elastic equations as a first-order approximation

Eq. (19)

σik,elastic=E1+νϵik+E1+νν12νϵjjδik[Po+EαΔT3(12ν)]δik=E1+νϵikPmδikPm=13σjj=E1+νν12νϵjj[Po+EαΔT3(12ν)]+E3(1+ν)ϵjj,
where E is the elastic modulus of the biotissue; ν is the Poisson’s ratio; ϵik is the strain tensor, which is a function of the displacement field u as (ui,j+uj,i)/2; Po is the environmental pressure, which is 1 atm in this study; and Pm is the mechanical pressure. Since the most water-rich biotissues (e.g., the cornea) are highly incompressible with ν close to 0.5, the mechanical pressure Pm can be approximated with thermodynamics pressure P, which can be estimated with Eq. (2) as

Eq. (20)

Pm=E1+νν12νϵjj[Po+EαΔT3(12ν)]+E3(1+ν)ϵjjsmall compared withother terms when  ν0.5E1+νν12νϵjj[Po+EαΔT3(12ν)]=P.

To ensure a numerically stable simulation, artificial viscous force is added in the simulation to minimize the artificial noises, which can be significant in this normal pressure driven hydrodynamics involving large pressure gradient and density variations.21 Adjusting the artificial viscous force to be one to two orders less than the pressure force to prevent alterations in the physical dynamics of the microcavitation. As a result, Eq. (21) includes the adopted stress tensor of biotissue in elastic solid phase

Eq. (21)

σik,elasticPδik+E1+νϵik+μart(Vi,k+Vk,i)artificial viscous term,
where μart is the artificial viscosity.

When the biotissue is in liquid or vapor phases, assumption of Newtonian fluid with ·V0 (when the material highly incompressible) is applied with the following stress tensor expression

Eq. (22)

σik,fluidPδik+μ(Vi,k+Vk,i)+μart(Vi,k+Vk,i)artificial viscous term,
where μ is the dynamic viscosity of the biotissue. When the biotissue is in a liquid or vapor phase, an artificial viscous force term is included in order to improve the numerical stability.

2.5.

Energy Conservation for the Water-Rich Biotissue

The local energy conversion in the biotissue is expressed as

Eq. (23)

et+(V·)e·(κT)+Qr,
where κ is the thermal conductivity of the biotissue, and Qr refers to the heat generation rate due to laser light absorption and the associated Joule heating. The Qr is determined from wave optics by solving a flat-top laser light propagating through a lens and then focused at a specified location inside the biotissue with time harmonic Maxwell’s equations as

Eq. (24)

·D=0;·B=0;×E=jωμoH;×H=jωϵrϵoE,
with ω is the laser frequency; E, D, B, and H are the electric field, electric displacement, magnetic induction, and field in the frequency domains, respectively; ϵr is the relative permittivity of the biotissue [note that absorption coefficient can be expressed as a function of ϵr as 4π×IMAG(ϵr)/λ with IMAG an operator determining the imaginary part of a complex number, and λ is the wavelength]; ϵo is the permittivity of free space; and μo is the permeability of the biotissue. The induced current density in the frequency domain is J=jωD. The time average joule heating ϕ per unit volume is ϕ=1/2Re(jωE·D*). The heat generation rate Qr is equal to ϕ during the laser pulse and is zero after the laser pulse. Figure 2(b) shows the simulation domain and the obtained spatial distribution of joule heating ϕ under time harmonic condition. The incident flat-top beam had a diameter of 200  μm with a wavelength (λ) of 4  μm. The hemispherical lens had a diameter of 300  μm with a refractive index of 2.5. The numerical aperture of the focal beam was 0.9. The focal position was 80  μm beneath the surface of the tissue, which is a common location during a LASIK surgery. The focal spot was 4  μm in diameter. The refractive index of the tissue domain is 1.37 with an absorption coefficient of 8759  m1. Using a long wavelength with a large absorption coefficient provides thermal microcavitation without avalanche breakdown in the biotissue. Furthermore, a body of resolution was used in order to reduce the original three-dimensional simulation to a quasi-two-dimensional (2-D) simulation.

Fig. 2

(a) Simulation domain for the plane wave propagation through a hemispherical lens and then focuses inside the biotissue. (b) Spatial distribution of joule heating inside the biotissue. A rainbow color scale is applied to illustrate the relative intensity value. The domain size of the joule-heating plot is 100  μm×150  μm.

JBO_22_5_056002_f002.png

3.

Results and Discussion

3.1.

Verification of the Integrated Analysis

To verify the capability of the developed integrated analysis during ns laser-induced thermal microcavitation in a water-rich biotissue, we first examined a 10-ns λ=4  μm laser pulse focused at 80  μm beneath the surface of the tissue with elastic modulus 60  kPa and Poisson ratio=0.4995 (i.e., highly incompressible tissue, which is true for water-rich tissues). The distribution of laser-induced joule heating inside the tissue is the same as that in Fig. 2(b). A 0.02-mJ laser energy was delivered from the incident plane of Fig. 2(a), in order to raise the tissue temperature and create a strong thermal microcavitation event. Since the dynamic process inside the biotissue is axi-symmetric, a quasi-2-D simulation was adopted to solve the microcavitation processes. The width and height of the quasi-2-D simulation domain were 100 and 200  μm, respectively. The initial temperature of the biotissue was 293.15 K, and the outer boundaries of the simulation domain were fixed at 293.15 K during the simulation. The top surface boundary was not fixed and maintained at 1 atm. The movement of the top surface was determined from the velocity of the boundary determined with an arbitrary Lagrangian Euler (ALE) scheme.22 ALE is a finite-element formulation, in which the mesh system can move arbitrary to trace the distorting boundary. The obtained results based on equation of state listed in Sec. 2.1, expression of internal energy listed in Sec. 2.2, and governing equations listed in Secs. 2.32.5 are found in Fig. 3.

Fig. 3

(a) Temperature increment (on top of 293.15 K with unit of Kelvin) inside the bio tissue at different delay times from the beginning of the 10-ns laser pulse. (b) Pressure distribution (with unit of Pa) inside the biotissue at different delay times. (c) Vapor quality inside the biotissue at different delay times. “1” means pure vapor phase.

JBO_22_5_056002_f003.png

As illustrated in the temperature plots of Fig. 3(a), the tissue temperature increases monotonically during the 10-ns laser pulse. However, the temperature increment is not exactly linear with respect to time due to the varying specific heat as a function of temperature and specific volume. The tissue cools down after the laser pulse with a speed much slower than the heating. The heat-affected zone also increases monotonically with respect to the delay time. The radius of the heat-affected zone is 25  μm at 1  μs from the beginning of the laser pulse. High pressure occurs around the focal spot during the laser pulse due to the rapid heating of the biotissue as shown in Fig. 3(b). The maximum pressure is 5×108  Pa at the end of the laser pulse, which is within the failure limit of most water-rich biotissue (>GPa). The thermally induced high pressure around the laser pulse causes a pressure wave propagating from the laser spot to the surrounding area seen in the pressure plots from 10 to 100  ns from the beginning of the laser pulse. The pressure front is nearly spherical. The pressure wave leaves the right boundary at t60  ns, which correspond to an average speed >1600(m/s). This propagation speed is much higher than the elastic wave speed in the biotissue, which is E(1ν)/ρ/(1+ν)/(12ν)282.75(m/s) when the tissue density is 1000(kg/m3). Note that no shear wave was observed in the simulation, which indicates that the microcavitation dynamics are dominated by normal pressure. Note that compared with the temperature distribution, the pressure shows a uniform distribution at the heat affected zone especially when the time is going on. The uniform distribution is attributed to the pressure waves, which can even out the pressure gradient in the heat-affected zone.

Figure 3(c) shows the evolution of the vapor quality, which is defined as the mass fraction of vapor inside the tissue during the microcavitation process. It is interesting to point out that the vapor bubble does not appear right after the end of the laser pulse when the maximum temperature occurs. Instead, the vapor bubble [i.e., red region in Fig. 3(c)] was observed at 16  ns from the beginning of the laser pulse, which is 6 ns after the end of the 10-ns laser pulse. This delay relates to the pressure wave propagation: the high temperature at the laser spot during the laser pulse does not convert the portion of tissue into a vapor bubble directly due to the accompanied high pressure around the laser spot. To trigger the vaporization process, an outgoing pressure wave is required to release a portion of the high pressure around the laser pulse to the surrounding, and it takes a few ns for the pressure wave to detach from the highly heated zone around the laser spot. The shape of the vapor bubble follows the joule heating profile [Fig. 2(b)] initially and rapidly evolves into a larger more spherical bubble within 1  μs. The averaged expansion speed of the vapor bubble (defined as the radius increasing rate in this study) is <20(m/s), which is much less than the elastic wave speed in the biotissue. Another interesting phenomenon observed in the vapor quality plot is the thickness of the phase transition zone (where the quality number varies from zero to one), which is very thin compared with the diameter of the vapor bubble once the diameter of the vapor bubble is larger than a couple of microns. This observation indicates that previous microcavitation simulations treating vapor and condensed phases as two separated domains with different governing equations without considering their transition is valid as long the vapor bubble is more than a few microns in diameter.

Based on the obtained temperature, pressure, and vapor quality plots, the integrated analysis reproduces the rapid heating, phase change, bubble formation, and pressure wave propagation at the beginning of a microcavitation process, which are the main characteristics observed in previous experiments. This qualitative consistency verifies the capability of the new integrated analysis in directly simulating a thermal microcavitation process in a water-rich biotissue when avalanche breakdown is not important and no significant plasma forms.

Figure 4 shows the spatial distribution of the specific heat cv of the biotissue obtained with equations of Sec. 2.2 at 50 ns from the beginning of the laser pulse. The cv value of the tissue is 4000  J/kg/K when it is far away from the laser spot, and the tissue remains in elastic solid phase. The cv value changes to 2000  J/kg/K around the laser spot when the biotissue is fully vaporized, and the specific heat is close to that of water vapor. Between these two pure phases, cv changes vigorously in the phase change zone and is not only a function of temperature but also a strong function of specific volume. The maximum cv (>4000  J/kg/K) happens when the vaporization of the water-rich biotissue is about to finish. Without considering the rapid change of cv, the simulation can be numerically unstable at locations where the phase change is about to start or stop; especially, when the temperature of the water-rich biotissue is close to the critical temperature of the water (647  K).

Fig. 4

Spatial distribution of specific heat cv of the biotissue with unit of J/kg/K at 50 ns from the beginning of the laser pulse. Distribution of the cv value at different delay times shows similar characteristics.

JBO_22_5_056002_f004.png

3.2.

Contribution of Laser Energy

After the verification of the capability of the developed integrated analysis, the analysis is conducted under the same condition as Sec. 3.1, except the laser energy is varied to 0.005, 0.01, and 0.015 mJ, respectively, to study the contribution of laser energy on the specific type of thermal macrocavitation process in a water-rich biotissue. Figure 5 shows these results.

Fig. 5

Temperature increment (on top of 293.15 K with unit of Kelvin), pressure distribution, and the vapor quality inside the water-rich biotissue at different delay times from the beginning of the 10-ns laser pulse. The laser energy labeled on the top of each group of figures is the energy delivered from the incident plane of Fig. 2(a).

JBO_22_5_056002_f005.png

As shown in temperature plots in Figs. 3 and 5, the peak temperature of the tissue appearing at the end of the laser pulse (i.e., 10 ns in this section) increases monotonically with respect to the incident laser energy. However, this peak temperature is more than doubled when the laser energy is doubled, which can be attributed to the variation of the specific heat of the tissue as a function of temperature and specific volume as described in Sec. 3.1. Another observation from the thermal plots is the variation of size and shape of the heat affected zone at each specific delay time when the laser energy is changed. For lower laser energies, little to no vaporization occurs in the water-rich biotissue, even at longer times after the laser pulse. As a result, the shape of the heat-affected zone stays the same as the shape of the laser-induced joule heating [c.f., Fig. 2(b)] at each delay time, which is not spherically symmetric. When the laser energy is more than the microcavitation threshold (e.g., 0.015 to 0.02 mJ for the laser and tissue condition studied in Secs. 3.1 and 3.2), the water-rich biotissue is significantly vaporized around the focal spot after 10s of ns, as shown in the vapor quality plots. As a result, the heat-affected zone gradually evolves into a more spherically symmetric shape at longer delay time. This can explain the spherically shaped vapor bubble and heat-affected zone reported in the most previous microcavitation experiments when high laser energy is applied.

Furthermore, based on the pressure plots especially at 50 ns from the beginning of the laser pulse, the shockwave speed increases when the incident laser energy is increased. This trend qualitatively fits with the blast wave theory, which describes the blast wave speed as a function of deposited laser energy.23 Peak pressure of the biotissue happens at the end of the laser pulse (i.e., 10 ns in this case), which increases with respect to the laser energy. This high pressure releases rapidly by the shockwave propagation. Note that due to the swelling of the biotissue at the top surface of the simulation domain, which prevents a significant compression of the surrounding biotissue during the expansion of the vapor bubble, the pressure of the surrounding biotissue maintains at 1  atm even at 1  μs from the beginning of the laser pulse.

3.3.

Contribution of the Elastic Modulus of the Biotissue

It was observed by Vogel that the elastic modulus of the biotissue can affect the oscillation period of the vapor bubble from the microcavitation process long time (>10  s of μs) after the laser pulse.1,24 Therefore, we expect that the elastic modulus may also affect the formation process of the microcavitation during and shortly after the laser pulse. To verify this hypothesis, laser and tissue conditions of Sec. 3.1 were repeated in this Sec. 3.3 simulation except the elastic modulus is increased from 60 to 120 kPa and 240 kPa, respectively. The corresponding results are listed in Fig. 6.

Fig. 6

Temperature increment (on top of 293.15 K with unit of Kelvin), pressure distribution, and the vapor quality inside the water-rich biotissue at different delay times from the beginning of the 10-ns laser pulse. Two different elastic modulus of the tissue, namely 120 and 240 kPa, are adopted in this portion of analysis.

JBO_22_5_056002_f006.png

When we compare the results in Figs. 6 and 3, it is interesting to see that the peak temperature and pressure at the end of the 10-ns laser pulse show minimum change when the elastic modulus is doubled or quadrupled. The propagation speed of the induced shockwave also shows little to no change when the elastic modulus is varied in this specific range. In addition, the sizes of the vapor bubble are almost the same for tissue with these three different elastic modulus. The insensitivity of the formation stage of the microcavitation processes to the range of the elastic modulus we tested can be understood from the driving mechanisms of each phenomenon during and shortly after the laser pulse. The peak temperature and pressure at the end of the laser pulse are determined by the specific heat and the equation of the state of the biotissue during the laser pulse. As long as the temperature of the tissue is higher than the dissociation temperature (340  K in this study), the specific heat and the equation of the state of the water-rich biotissue are not affected by the elastic modulus of the tissue. This explains the independence of the peak temperature and pressure with the elastic modulus of water-rich biotissue under the higher laser energy condition we simulated. The propagation speed of the shockwave (>1600  m/s), when it is much higher than the elastic wave speed of the water-rich biotissue (i.e., 282.75, 399.87, and 565.50  m/s for elastic modulus equal to 60, 120, and 240 kPa), can be described with blast wave theory. In blast wave theory, the propagation speed of the shockwave is only determined with the absorbed laser energy, density of the surrounding medium, and the delay time after the laser pulse.23 As a result, the propagation speed of the strong shockwave (i.e., the blast wave) during the first hundreds of ns showed little dependence on the elastic modulus of the water-rich biotissue.

On the other hand, the expansion of the vapor bubble during the beginning of the microcavitation process is due to the outgoing thermal energy diffusion from the laser spot and the resulting phase change of the biotissue. When biotissue changes into vapor phase after receiving enough thermal energy, the associated pressure increment can further expand the size of the vapor bubble to the surrounding tissue. Since the change of the elastic modulus of the biotissue does not affect the peak temperature and the thermal diffusion coefficient of the biotissue, the expansion speed of the vapor bubble does not show a significant change during the first 1  μs as illustrated in the result of the simulation when the elastic modulus is doubled or quadrupled. Furthermore, we assume that the top boundary free deformation of the simulation domain in our simulation prevents an accumulation of the pressure field in the surrounding biotissue during the bubble expansion. Furthermore, this helps to minimize the contribution of the elastic modulus on the expansion process of the vapor bubble during the early stage of the microcavitation process.

3.4.

Contribution of Laser Pulse Duration

Based on the understanding that the variation of the specific heat during the laser–tissue interaction can significantly affect the microcavitation process (as observed in Sec. 3.2), we expect the laser pulse duration to also affect the microcavitation process. When the laser pulse duration is longer, the pressure around the laser spot due to the temperature rise has time to release during the laser pulse, and the tissue can be converted (or partially converted) into a vapor phase, which has lower specific heat compared with the condensed phase. As a result, we expect the peak temperature of the biotissue at the end of the laser pulse to increase with longer laser pulses under the same laser energy. Figure 7 shows the simulation result for the same laser and tissue condition as case 3a, except the pulse duration increased from 10 to 20 ns. As we expected, the peak temperature increment of the tissue raised from 1373 to 1759 K at the end of the laser pulse when the pulse duration increased from 10 to 20 ns while the total input laser energy remains the same. The higher temperature also causes a higher pressure inside the biotissue. Owning to the higher temperature and pressure, the vapor bubble experiences a higher expansion speed, which results in a larger vapor bubble at 1  μs from the beginning of the laser pulse when the pulse duration is increased from 10 to 20 ns. Though the bubble size is increased, due to the much lowered pressure increase rate when the pulse duration is doubled, the induced shockwave has a lower speed when the pulse duration is increased from 10 to 20 ns.18

Fig. 7

Temperature increment (on top of 293.15 K with unit of Kelvin), pressure distribution, and the vapor quality inside the water-rich biotissue at different delay times from the beginning of the 20-ns laser pulse.

JBO_22_5_056002_f007.png

4.

Conclusion

Even with numerous experimental studies, theoretical analysis describing the formation stage of a microcavitation process in a biotissue is still in a premature stage. The rapid phase change and the associated physical property variations cause obstacles in simulating the microcavitation process, especially during its formation stage. In this study, we addressed these issues by presenting the pressure field of the biotissue in different phases as a function of local pressure and density with an equation of state. As a first approximation, we assumed the equation of state of a water-rich biotissue can be approximated with the equation of state of water when it is in the liquid and vapor phases. When the biotissue is in the elastic material phase, a small modification with a scaling factor is introduced to the equation of state of water to satisfy experimentally measured bulk modulus (i.e., inverse of the compressibility) of the biotissue. With the analytical expression of the pressure field, we also derived an analytical expression of the specific heat of the biotissue as a function of the local temperature and specific volume. After integrating the analytical expressions of the pressure and specific heat of the biotissue in different phases with mass, momentum, and energy conservation equations, we successfully reproduced the laser-induced heating. We also successfully reproduced the vapor bubble formation, expansion, and shockwave propagation inside a water-rich biotissue during and shortly after the laser pulse during a thermal microcavitation process with no avalanche breakdown and minimal plasma formation. Such a specific type of thermal microcavitation process is common in laser surgeries when a long wavelength nanosecond pulsed laser (e.g., Er-YAG ns laser) is utilized.

With this new analysis, we identified several characteristics of the thermal microcavitation process in a water-rich biotissue during and shortly after a laser pulse. First, the temperature and pressure fields of biotissue show a nonlinear change with respect to the incident laser energy during the microcavitation process. This nonlinear change is mainly due to the variation of the specific heat as a function of local temperature and specific volume. Second, the change of laser pulse duration can also affect the temperature and pressure fields of the biotissue during the specific thermal microcavitation process. Finally, the temperature, pressure, and propagation speed of laser-induced shockwaves in a water-rich biotissue show little dependence on the elastic modulus of the biotissue during the specific microcavitation process, especially when it is close to an open boundary.

Disclosures

No conflicts of interest, financial or otherwise, are declared by the authors.

Acknowledgments

This research was supported by the Air Force Contract FA8650-14-D-6519 (Engility Corporation) and the US Air Force Summer Faculty Fellowship Program.

References

1. 

E. A. Brujan and A. Vogel, “Stress wave emission and cavitation bubble dynamics by nanosecond optical breakdown in a tissue phantom,” J. Fluid Mech., 558 281 –308 (2006). http://dx.doi.org/10.1017/S0022112006000115 Google Scholar

2. 

J. Noack and A. Vogel, “Laser-induced plasma formation in water at nanosecond to femtosecond time scales: calculation of thresholds, absorption coefficients, and energy density,” IEEE J. Quantum Electron., 35 1156 –1167 (1999). http://dx.doi.org/10.1109/3.777215 Google Scholar

3. 

A. Vogel et al., “Femtosecond-laser-induced nanocavitation in water: implications for optical breakdown threshold and cell surgery,” Phys. Rev. Lett., 100 038102 (2008). http://dx.doi.org/10.1103/PhysRevLett.100.038102 Google Scholar

4. 

N. Tinne et al., “Interaction dynamics of spatially separated cavitation bubbles in water,” J. Biomed. Opt., 15 068003 (2010). http://dx.doi.org/10.1117/1.3526366 Google Scholar

5. 

J. R. Blake and D. C. Gibson, “Cavitation bubbles near boundaries,” Annu. Rev. Fluid Mech., 19 99 –123 (1987). http://dx.doi.org/10.1146/annurev.fl.19.010187.000531 Google Scholar

6. 

T. Kodama and Y. Tomita, “Cavitation bubble behavior and bubble-shock wave interaction near a gelatin surface as a study of in vivo bubble dynamics,” Appl. Phys. B: Lasers Opt., 70 139 –149 (2000). http://dx.doi.org/10.1007/s003400050022 Google Scholar

7. 

M. S. Schmidt et al., “Temperature dependence of nanosecond laser pulse thresholds of melanosome and microsphere microcavitation,” J. Biomed. Opt., 21 015013 (2016). http://dx.doi.org/10.1117/1.JBO.21.1.015013 Google Scholar

8. 

M. S. Schmidt et al., “Temperature dependence of melanosome microcavitation thresholds produced by single nanosecond laser pulses,” Proc. SPIE, 9321 932108 (2015). http://dx.doi.org/10.1117/12.2079934 Google Scholar

9. 

M. S. Schmidt et al., “Trends in nanosecond melanosome microcavitation up to 1540 nm,” J. Biomed. Opt., 20 095011 (2015). http://dx.doi.org/10.1117/1.JBO.20.9.095011 Google Scholar

10. 

C. R. Thompson et al., “Melanin granule model for laser-induced thermal damage in the retina,” Bull. Math. Biol., 58 809 –810 (1996). http://dx.doi.org/10.1007/BF02459483 BMTBAP 0092-8240 Google Scholar

11. 

B. S. Gerstman et al., “Laser induced bubble formation in the retina,” Laser Surg. Med., 18 10 –21 (1996). http://dx.doi.org/10.1002/(ISSN)1096-9101 Google Scholar

12. 

M. Jaunich et al., “Bio-heat transfer analysis during short pulse laser irradiation of tissues,” Int. J. Heat Mass Transfer, 51 5511 –5521 (2008). http://dx.doi.org/10.1016/j.ijheatmasstransfer.2008.04.033 Google Scholar

13. 

R. P. Godwin et al., “Aspherical bubble dynamics and oscillation times,” Proc. SPIE, 3601 225 (1999). http://dx.doi.org/10.1117/12.350042 PSISDG 0277-786X Google Scholar

14. 

E. Faraggi, B. S. Gerstman and J. M. Sun, “Biophysical effects of pulsed lasers in the retina and other tissues containing strongly absorbing particles: shockwave and explosive bubble generation,” J. Biomed. Opt., 10 064029 (2005). http://dx.doi.org/10.1117/1.2139970 Google Scholar

15. 

J. M. Sun, B. S. Gerstman and B. Li, “Bubble dynamics and shock waves generated by laser absorption of a photoacoustic sphere,” J. Appl. Phys., 88 2352 –2362 (2000). http://dx.doi.org/10.1063/1.1288507 Google Scholar

16. 

B. S. Gerstman, “Nanosecond laser pulses: from retinal damage to refined surgical treatment of congenital nevi and melanoma,” Proc. SPIE, 2975 180 (1997). http://dx.doi.org/10.1117/12.275474 PSISDG 0277-786X Google Scholar

17. 

W. B. Telfair et al., “Histological comparison of corneal ablation with Er: YAG laser, Nd: YAG optical parametric oscillator, and excimer laser,” J. Refractive Surg., 16 40 –50 (2000). Google Scholar

18. 

V. P. Carey, Statistical Thermodynamics and Microscale Thermophysics, Cambridge University Press, Cambridge, United Kingdom (1999). Google Scholar

19. 

T. H. Ahmed, Equations of State and PVT Analysis: Applications for Improved Reservoir Modeling, 2nd ed.Gulf Professional Publishing, an imprint of Elsevier, Cambridge, Massachusetts (2016). Google Scholar

20. 

W. Wagner, H.-J. KretzschmarW. Wagner, International Steam Tables: Properties of Water and Steam Based on the Industrial Formulation IAPWS-IF97: Tables, Algorithms, Diagrams, and CD-ROM Electronic Steam Tables: All of the Equations of IAPWS-IF97 Including a Complete Set of Supplementary Backward Equations for Fast Calculations of Heat cycles, Boilers, and Steam Turbines, 2nd ed.Springer, Berlin (2008). Google Scholar

21. 

D. A. Anderson, R. H. Pletcher and J. C. Tannehill, Computational Fluid Mechanics and Heat Transfer, 3rd ed.CRC Press, Taylor & Francis Group, Boca Raton (2013). Google Scholar

22. 

E. Stein, R. De Borst and T. J. R. Hughes, Encyclopedia of Computational Mechanics, John Wiley, Chichester, West Sussex (2004). Google Scholar

23. 

L. I. Sedov, Similarity and Dimensional Methods in Mechanics, 10th ed.CRC Press, Boca Raton, Florida (1993). Google Scholar

24. 

E. A. Brujan et al., “Dynamics of laser-induced cavitation bubbles near an elastic boundary,” J. Fluid Mech., 433 251 –281 (2001). http://dx.doi.org/10.1017/S0022112000003347 Google Scholar

Biographies for the authors are not available.

© 2017 Society of Photo-Optical Instrumentation Engineers (SPIE) 1083-3668/2017/$25.00 © 2017 SPIE
Sy-Bor Wen, Kevin J. Ly, Arun Bhaskar, Morgan S. Schmidt, and Robert J. Thomas "Direct numerical simulation of the initial stage of a thermally induced microcavitation in a water-rich biotissue triggered by a nanosecond pulsed laser," Journal of Biomedical Optics 22(5), 056002 (3 May 2017). https://doi.org/10.1117/1.JBO.22.5.056002
Received: 24 January 2017; Accepted: 11 April 2017; Published: 3 May 2017
Lens.org Logo
CITATIONS
Cited by 4 scholarly publications and 1 patent.
Advertisement
Advertisement
RIGHTS & PERMISSIONS
Get copyright permission  Get copyright permission on Copyright Marketplace
KEYWORDS
Pulsed laser operation

Laser tissue interaction

Laser energy

Laser vision correction

Tissues

Biological research

Liquids

Back to Top