COMSOL 6.4 - Focused Ultrasound Induced Heating in Tissue Phantom

Focused Ultrasound Induced Heating in Tissue Phantom

When an ultrasound beam passes through a volume of tissue, some of the energy of the primary acoustic field is absorbed locally by the tissue and turned into heat. This results in a temperature increase whose magnitude is a function of the physical properties of the medium (acoustic absorption coefficient, density, and specific heat), the properties of the ultrasound device (beam geometry), and the frequency and time-averaged acoustic intensity of the acoustic field; see Ref. 1 and Ref. 2. The actual temperature rise that can be obtained also depend on the conduction and convection properties of the tissue involved, for example, the blood perfusion rate of the tissue.

Figure 1: Focused ultrasound makes selective and targeted heating possible: heating of tissues lying within the focal volume can be achieved with minimal damage to nearby healthy tissue and other structures lying elsewhere in the path of the beam.

This application requires the Acoustics Module and Heat Transfer Module.

Therapeutic applications of ultrasound typically involve focused beams, which allow the ultrasound energy to be directed into a small area within the tissue region that needs the treatment. Moreover, heating of tissues lying within the focal volume can be achieved with minimal damage to nearby healthy tissue and other structures lying elsewhere in the path of the beam, as shown in Figure 1.

Depending on the dosage parameters — that is, field intensity and exposure time — the clinical applications of focused ultrasound (FU) are generally grouped into two categories, namely “ultrasound hyperthermia” and “focused ultrasound surgery” (FUS); see Ref. 3. Generally, during hyperthermia applications, tissues are exposed to ultrasound for long periods (from 10 to 60 minutes) at lower intensity levels, such that the irradiated tissue temperature is elevated and maintained at 41 °C to 45 °C during the therapy. The biological change thus induced can be reversed. In contrast, focused ultrasound surgery utilizes intense, relatively short bursts (0.1 s to 30 s) to induce irreversible changes in the focal tissue volume. In this type of applications, nonlinear acoustic effects and acoustic cavitation usually play essential roles; the tissue temperature in the focal zone can reach 70 °C to 90 °C within a few seconds. This technique is also known as ultrasound ablation.

The thermal effect of focused ultrasound also brings concerns about possible harmful effects of diagnostic ultrasound, especially in obstetrical examinations when the fetus is exposed to ultrasound. The U.S. Food and Drug Administration has set up regulations on the maximum thermal index, a dosage parameter reflecting the combined effect of temperature and exposure time that is an estimate of risk from heat (Ref. 4). Diagnostic ultrasound systems now come with displays with the Thermal Index (TI) and the Mechanical Index (MI), the other estimate of risk from the nonthermal effects of ultrasound in order to meet the U.S. government’s regulations. Safety in the use of ultrasound has also been addressed extensively in the academic field; see, for example, Ref. 5, Ref. 6, and Ref. 7.

The current model is inspired by the experiments to measure focused ultrasound induced heating in a tissue phantom from Ref. 8. The model uses the same geometry and material properties as in Ref. 8. The model exemplifies how to use COMSOL Multiphysics to model tissue heating induced by focused ultrasound. The simulation results are compared to the experimental data in the reference.

Model Definition

Figure 2 shows the geometry simulated in this model. Both the tissue phantom and the acoustic transducer are immersed in water. The transducer is bowl shaped with a focal length of 62.64 mm, an aperture of 35 mm in radius, and a hole of 10 mm in radius in the center. The tissue phantom has the shape of a cylinder with 53.6 mm in radius and 80.5 mm in length. The tissue phantom and the transducer are arranged coaxially so the model can be defined as being 2D axisymmetric.

The transducer is driven at the frequency of 1 MHz. It is turned on for 1 second and then turned off to let the tissue phantom cool down completely in water. The model thus solves for the heating of the tissue phantom for 1 second and then simulates the cooling process after the acoustic source is turned off.

Figure 2: The geometry of the acoustic transducer and tissue phantom. The transducer is bowl shaped with a hole in the center. Both the tissue phantom and the transducer are immersed in water. The axisymmetric geometry allows for a 2D axisymmetric simulation.

The model uses the Pressure Acoustics, Frequency Domain interface to model the stationary acoustic field in the water and the tissue domain, to obtain the acoustic intensity distribution in the tissue phantom. The absorbed acoustic energy is calculated and used as the heat source for the Bioheat Transfer interface model. Because the acoustic focal region (like the heated area) is much smaller than the size of the tissue phantom, the thermal simulation is performed only in the tissue domain.

The wave equation solved is the homogeneous Helmholtz equation in 2D axisymmetric cylindrical coordinates:

(1)

Here r and z are the radial and axial coordinates, p is the acoustic pressure, and ω is the angular frequency. The density, ρc, and the speed of sound, cc, are complex-valued to account for the material’s damping properties.

Using Equation 1 involves the assumption that the acoustic wave propagation is linear and also that the amplitude of shear waves in the tissue domain are much smaller than that of the pressure waves. Nonlinear effects and shear waves are therefore neglected.

Given the acoustic pressure field, the acoustic intensity field is readily derived. The heat source Q for thermal simulation, given in the plane-wave limit, is then calculated as

(2)

where αABS is the acoustic absorption coefficient, I is the acoustic intensity magnitude, p is the acoustic pressure, and v is the acoustic particle velocity vector. In COMSOL, the intensity is a derived variable where the magnitude can be accessed as acpr.I_mag. Likewise, the dissipated power density (for plane waves) is defined as acpr.Q_pw. The heat source Q is thus readily calculated once the acoustic field is solved.

For further details about the intensity variables choose > and then search for the string Modeling with the Pressure Acoustics Branch (FEM-Based Interfaces). Then select the section Postprocessing Variables. Two parts exist here; one describing the intensity variables and one describing power dissipation variables.

Inserting the volumetric acoustic heat source into the Pennes’ Bioheat Transfer equation to model heat transfer within biological tissue gives

where T is the temperature, ρ is the density, Cp is the specific heat, k is the thermal conductivity, ρb is the density of blood, Cb is the specific heat of blood, wb is the blood perfusion rate, Tb is the temperature of the blood, Q is the heat source (the absorbed ultrasound energy calculated from Equation 2), and Qmet is the metabolic heat source.

In this model, assume that the tissue properties do not change when the temperature rises. Blood perfusion is also neglected.

Figure 3 shows the model geometry and material domains defined in the model. Four cylindrical perfectly matched layers (PMLs) (r1-r4) and one spherical PML (c1) are used to absorb the outgoing waves. The pressure acoustics simulation is performed in all domains while the heat transfer model is only applied in the tissue phantom domain.

Figure 3: Model geometry. Water domains are shown in blue and tissue phantom domains are shown in purple. Four cylindrical PMLs (r1-r4) and one spherical PML (c1) are used to absorb the outgoing waves.

To accurately resolve the sharp pressure gradient in the focal region, the model uses a fine mesh with size λ/8 (where λ is the wavelength) within that region (both for the heat and the acoustics problem). A coarser mesh with size λ/5 in water is used for the other acoustic domains. Both the acoustic pressure and the temperature uses the default quadratic (2nd order) elements as discretization.

Table 1 shows the material properties used in the model simulation. The properties used for the tissue phantom are the measured data described in Ref. 8. For comparison, the table also lists properties for human tissue published in Ref. 9.

Table 1: Material properties used in the model.
Property	Density (kg/m3)	Speed of sound (m/s)	Attenuation (Np/m/MHz)	Specific heat (J/(kg·K))	Thermal conductivity (W/(m·K))
Water (at 293.7 K)	1000	1483	0.025	N/A	N/A
Tissue phantom	1044	1568	8.55	3710	0.59
Human tissue	1000–1100	1450–1640	4.03–17.27	3600–3890	0.45–0.56

Results and Discussion

Figure 4 depicts the acoustic pressure. The ultrasonic beam travels through the water layer and into the tissue phantom. The beam converges into a focal area where the acoustic pressure amplitude reaches as high as 1.11 MPa at the focal point. This result agrees well with the results presented in Ref. 8. Note the diffraction-like pattern near the edges (especially at the top of the computational domain); this is not a physical phenomenon but an effect of the perfectly matched layer.

Figure 4: The acoustic pressure field in the water and tissue domains.

A depiction of the acoustic intensity magnitude is given in Figure 5. This plot shows more clearly how the acoustic energy is focused and distributed in the area of interest. Most of the heating happens in the oval-shaped focal area which is about 8 mm long and 1.3 mm wide. Figure 6 shows the acoustic pressure amplitude profile along the z-axis (r = 0). By zooming in around the peak pressure amplitude in Figure 6, the exact location of the acoustic focus is seen to be on-axis and 35 mm away from the tissue phantom and water interface (at z = 59.6 mm). Figure 7 shows the acoustic pressure amplitude profile along the radial direction in the focal plane (z = 59.6 mm).

Figure 8 shows the result of the temperature rise in the tissue phantom for a focal pressure amplitude of 1.11 MPa after 1 second of insonation. At t = 1 s, the maximum temperature rise is about one degree. The oval-shaped heated spot is about the same size as that of the acoustic focal area, which is more easily seen in the contour plot of the temperature as shown in Figure 9.

Figure 5: Surface plot of the acoustic intensity field, showing the acoustic energy to be concentrated to the focal region.

Figure 10 plots the heating and cooling curves at acoustic focus and 0.5 mm off acoustic focus in the focal plane. This is again for a focal pressure amplitude of 1.11 MPa. At these locations, the tissue heats up when insonated and then cools down through natural conduction. The result at the acoustic focus agrees well with the results presented in Ref. 8.

This example shows how to model tissue heating induced by focused ultrasound when the acoustic pressure at focus is well below acoustic cavitation threshold. Because the model does not take nonlinearity into account it can be solved in the frequency domain. Linear acoustic is valid in the limit where p << ρc2. The maximal pressure in the focal region is much smaller than the product ρc2 = 2.6 109 Pa and thus the linear assumption of the frequency domain is valid. At pressures comparable to and above the limit, acoustic energy is pumped out from the fundamental frequency to higher harmonics (Ref. 10). This nonlinear behavior can be modeled in the time domain using the Nonlinear Acoustics (Westervelt) domain feature of the Pressure Acoustics, Transient interface. In tissue, the absorption coefficient is close to f 1.1 (for water, the power law for absorption is f 2). Higher harmonics result in a narrower focal region and more energy dissipation. Therefore at higher sound pressure levels, the current model would tend to underestimate the total energy deposition associated with the absorption of the ultrasonic wave and consequently also the maximum temperature rise at the focal region.

Figure 6: Acoustic pressure amplitude profile along the symmetry axis (r = 0).

Figure 7: Acoustic pressure amplitude profile along the radial direction in the focal plane.

Figure 8: Surface plot of the temperature rise in the tissue phantom after 1 second insonation for a focal pressure amplitude of 1.11 MPa.

Figure 9: Contour plot of the temperature rise in the tissue phantom after 1 second of insonation for a focal pressure amplitude of 1.11 MPa. The oval-shaped heated spot is about the same size as that of the acoustic focal area.

Figure 10: Heating and cooling curves at acoustic focus and 0.5 mm off-axis in the focal plane for 1 second of insonation with a focal pressure amplitude of 1.11 MPa.

The effect of energy dissipation serves to counteract the nonlinear effect and thus mitigate the waveform distortion (Ref. 10). Therefore, the assumption of linear progressive wave motion should remain good as long as the nonlinearity is relatively small, especially for hyperthermia applications. The simulation results show that this model provides a good estimate of both acoustic field and temperature rise for focal pressure up to 1.11 MPa. It not only simulates the heating and cooling behavior at the focal region but also gives other useful information, such as the shape of the heated area and the sidelobe heating effect at those surrounding locations with pressure maxima outside the focal region (which is the main lobe). The results also show that the sidelobes in the intensity field are mitigated in the temperature response due to the smoothing effect of conduction, as shown in Figure 11. In general, this model suggests that the temperature change is roughly proportional to the acoustic intensity.

Figure 11: The temperature and intensity profiles show that the sidelobes in the intensity field are mitigated in the temperature response due to the smoothing effect of the thermal conduction.

Notes About the COMSOL Implementation

This model uses two studies. Study 1 solves for Pressure Acoustics in the Frequency Domain to obtain the acoustic pressure field. In Study 2, the heat source calculated from the solution of Study 1 is then used in a time-dependent analysis to model heating and cooling of the tissue phantom.

References

1. K.J. Parker, “The Thermal Pulse Decay Technique for Measuring Ultrasonic Absorption Coefficients,” J. Acoust. Soc. Am., vol. 74, no. 5, pp. 1356–1361, 1983.

2. R.L. Clarke and G. ter Haar, “Temperature Rise Recorded during Lesion Formation by High Intensity Focused Ultrasound,” Ultrasound Med. Biol., vol. 23, no. 2, pp. 299–306, 1997.

3. N.T. Sanghvi, K. Hynynen, and F.L. Lizzi, “New Developments in Therapeutic Ultrasound,” IEEE Eng. Med. Biol. Mag., vol. 15, no. 6, pp. 83–92, 1996.

4. “Exposure Criteria for Medical Diagnostic Ultrasound: II. Criteria Based on All Known Mechanisms (Report No 140),” Diagnostic Ultrasound Safety: A Summary of the Technical Report, issued by the National Council on Radiation Protection and Measurements (NCRP), 2002.

5. P.N.T. Wells, “Biological Effects,” Biomedical Ultrasonics, chap. 9, Academic Press, New York, 1977.

6. G. ter Haar, “Biological Effects of Ultrasound in Clinical Applications,” Ultrasound: Its Chemical, Physical, and Biological Effects, K.S. Suslick, ed., VCH Publishers, New York, 1988.

7. W.L. Nyborg, “Biological Effects of Ultrasound: Development of Safety Guidelines. Part II: General Review,” Ultrasound Med. Biol., vol. 27, no. 3, pp. 301–333, 2001.

8. J. Huang, R.G. Holt, R.O. Cleveland, and R.A. Roy, “Experimental Validation of a Tractable Numerical Model for Focused Ultrasound Heating in Flow-through Tissue Phantoms,” J. Acoust. Soc. Am., vol. 116, no. 4, pt. 1, pp. 2451–2458, 2004.

9. F.A. Duck, Physical Properties of Tissue, Academic Press, 1990.

10. M.F. Hamilton and D.T. Blackstock eds., Nonlinear Acoustics, Academic Press, 1998.

Acoustics_Module/Ultrasound/ultrasound_induced_heating

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

In the tree, select > .

Click .

Click

In the tree, select > .

Click

Global Definitions

Define the parameters used in the model.

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
d0	3.8[nm]	3.8E-9 m	Displacement amplitude of transducer
z_tissue	24.6[mm]	0.0246 m	Starting position of tissue phantom
T0	293.7[K]	293.7 K	Initial temperature value
alpha_water	0.025[1/m]	0.025 1/m	Absorption coefficient of water
alpha_tissue	8.55[1/m]	8.55 1/m	Absorption coefficient of tissue phantom
f0	1[MHz]	1E6 Hz	Source frequency

Define a step function to turn off the acoustic source after 1 second of insonation. This is used for the transient thermal simulation. To improve convergence, define a smoothing transition zone that gently decreases the amplitude of the source to zero.

Step 1 (step1)

In the toolbar, click

and choose > .

In the window for , locate the section.

In the text field, type 1.

In the text field, type 0.

Click to expand the section. In the text field, type 0.005.

Geometry 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Now, proceed to create the geometry. A sketch is depicted in Figure 2 and Figure 3. First, create the bowl-shaped transducer used as the acoustic source.

Circle 1 (c1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 62.64.

Locate the section. In the text field, type 62.64.

Click

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 35.

In the text field, type 10.69.

Click

Intersection 1 (int1)

In the toolbar, click

and choose .

Click in the window and then press Ctrl+A to select both objects.

In the window for , click

Create the tissue phantom domain and add layers for the perfectly matched layer (PML).

Rectangle 2 (r2)

In the toolbar, click

In the window for , locate the section.

In the text field, type 48.6.

In the text field, type 75.5.

Locate the section. In the text field, type z_tissue.

Click to expand the section. In the table, enter the following settings:


Layer name	Thickness (mm)
Layer 1	5

Select the checkbox.

Clear the checkbox.

Click

Create the water domains and, also here, add layers for the PML.

Rectangle 3 (r3)

In the toolbar, click

In the window for , locate the section.

In the text field, type 48.6.

In the text field, type z_tissue-10.69.

Locate the section. In the text field, type 10.69.

Locate the section. In the table, enter the following settings:


Layer name	Thickness (mm)
Layer 1	5

Select the checkbox.

Clear the checkbox.

Click

Circle 2 (c2)

In the toolbar, click

In the window for , locate the section.

In the text field, type 15.

In the text field, type 90.

Locate the section. In the text field, type 0.80336.

Locate the section. In the text field, type -90.

Click to expand the section. In the table, enter the following settings:


Layer name	Thickness (mm)
Layer 1	5

Click

Union 1 (uni1)

In the toolbar, click

and choose .

Select the objects and only.

In the window for , click

Delete Entities 1 (del1)

In the window, right-click and choose .

On the object , select Boundary 10 only.

In the window for , click

Finally, create an oval-shaped focal region where a finer mesh can be applied to resolve high gradient in the pressure field.

Ellipse 1 (e1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 7.5.

In the text field, type 1.5.

In the text field, type 180.

Locate the section. In the text field, type z_tissue+35.

Locate the section. In the text field, type 270.

Click

Form Union (fin)

In the toolbar, click

Click the

button in the toolbar.

The geometry should look like the one depicted in the figure below.

In the window, click .

Definitions

Define the Perfectly Matched Layer domains for the Pressure Acoustics simulation. The rational stretching type is more efficient than the polynomial stretching type in these open problems, whereas the polynomial is preferred in waveguide-like problems.

Perfectly Matched Layer 1 (pml1)

In the window, expand the > node.

Right-click and choose .

Select Domains 7, 9, and 10 only.

It might be easier to select the domains by using the window. To open this window, in the toolbar click and choose . (If you are running the cross-platform desktop, you find in the main menu.)

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

Perfectly Matched Layer 2 (pml2)

Right-click and choose .

It is best to define a separate Perfectly Matched Layer domain for each different material present in the model.

In the window for , locate the section.

Click

Select Domain 8 only.

Perfectly Matched Layer 3 (pml3)

In the toolbar, click

Select Domain 1 only.

In the window for , locate the section.

From the list, choose .

Now, set up the physics for Pressure Acoustics, Frequency Domain.

Pressure Acoustics, Frequency Domain (acpr)

In the window, under click .

In the window for , locate the section.

From the list, choose .

Now, specify the physical properties for the water domains.

Pressure Acoustics 1

In the window, under > click .

In the window for , locate the section.

From the list, choose .

In the α text field, type alpha_water.

Add a second Pressure Acoustics Material Model node for the tissue domains. You will define the tissue material below.

Pressure Acoustics 2

In the toolbar, click

and choose .

Select Domains 5–7, 9, and 10 only.

In the window for , locate the section.

From the list, choose .

In the α text field, type alpha_tissue.

Define the normal displacement amplitude d0 at the surface of the ultrasound transducer.

Normal Displacement 1

In the toolbar, click

and choose .

Select Boundary 32 only.

In the window for , locate the section.

In the dn text field, type d0.

Now, specify the physics of the Bioheat Transfer interface.

Bioheat Transfer (ht)

In the window, under click .

In the window for , locate the section.

Click

Select Domains 5 and 6 only.

Initial Values 1

In the window, under > click .

In the window for , locate the section.

In the T text field, type T0.

Add the absorbed acoustic energy as the domain heat source. The variable acpr.Q_pw represents the dissipated power density for plane waves.

Heat Source 1

In the toolbar, click

and choose .

Select Domains 5 and 6 only.

In the window for , locate the section.

In the Q0 text field, type acpr.Q_pw*step1(t[1/s]-1).

Apply a constant temperature boundary condition on the computational boundaries of the thermal simulation.

Temperature 1

In the toolbar, click

and choose .

Select Boundaries 9, 14, and 20 only.

In the window for , locate the section.

In the T0 text field, type T0.

The next step is to set up the materials used in the model. Use the default water material and define your own tissue material.

Add Material

In the toolbar, click

to open the window.

Go to the window.

In the tree, select > .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Materials

Define the material properties of tissue phantom and apply it to the tissue phantom domains.

Tissue phantom

In the window, under right-click and choose .

In the window for , type Tissue phantom in the text field.

Select Domains 5–7, 9, and 10 only.

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	1044	kg/m³	Basic
Speed of sound	c	1568	m/s	Basic
Thermal conductivity	k_iso ; kii = k_iso, kij = 0	0.59	W/(m·K)	Basic
Heat capacity at constant pressure	Cp	3710	J/(kg·K)	Basic