High-Intensity Focused Ultrasound (HIFU) Propagation Through a Tissue Phantom

The high-intensity focused ultrasound is used in many different biomedical applications, such as thermal ablation of tumors, transcranial HIFU surgery, shock wave lithotripsy, etc. A HIFU signal is focused on a small focal zone, where its intensity reaches higher levels. In this case, nonlinear effects may become significant and one speaks of the propagation of finite amplitude sound waves. The nonlinear nature of the phenomenon results in the energy transfer from lower to higher harmonics. The contribution of the higher-order harmonics grows with the propagation distance.

This example shows how to model nonlinear propagation of HIFU through a dissipative media using the Nonlinear Pressure Acoustics, Time Explicit physics interface. The interfaces solves the system of nonlinear acoustic equations in the form of a hyperbolic conservation law, see Ref. 1, using the discontinuous Galerkin finite element method (dG-FEM) and an explicit time integration scheme. The used approach is suitable when the cumulative nonlinear effects, that is, those cumulated through propagation, dominate the local nonlinear effects.

In this model, the emitted signal is a tone burst pulse that occupies only a limited part of the computational domain on its way. Adaptive mesh refinement is used for automatic remeshing following the signal propagation. This ensures that the mesh is fine enough to resolve the higher-order harmonics where it is required.

Model Definition

The focusing of an ultrasonic signal is typically achieved either by using a phase delay or a focusing lens on the transducer side. In this model, a spherically focused ultrasound transducer with a concave lens is used to emit the signal. The transducer housing and the lens are assumed to be rigid. The model setup is shown in Figure 1. The model geometry is axially symmetric.

The spherical lens of radius r and aperture a emits a signal that is focused at the focal point F located in the tissue phantom. The signal is a five-cycle tone burst with the amplitude P0 = 0.1 MPa and the center frequency f0 = 1 MHz as shown in Figure 2. The signal amplitude at the source location is enough th generate of higher-order harmonics, but not high enough for the formation of shocks, and therefore no special shock-capturing techniques are required in this model.

Figure 1: Model geometry.

The traveling time of the signal to the focal point can be computed as

where c and d are the speed of sound and the traveling distance in the corresponding materials, respectively. At the frequency of 1 MHz, the time to focus is approximately equal to 40/f0.

Figure 2: Source signal.

Results and Discussion

The evolution of the signal traveling from the source is illustrated in Figure 3. The tone burst pulse has left the source at t = 10 μs and travels towards the boundary between the water and the tissue phantom domains, which it passes at t = 20 μs with some part being reflected back to the source. The focusing of the signal becomes visible at t = 30 μs and reaches its maximum at t = 40 μs.

Figure 3: Propagation of the ultrasonic signal at times t = 10, 20, 30, and 40 μs.

Figure 3 also shows that the signal strength increases as it comes closer to the focal zone. This is confirmed by the plots of the signal at the water/tissue interface and at the focal point depicted in Figure 4. The acoustic pressure amplitude at the focal point is about 10 times greater than that at the water/tissue interface. Another observation is that the positive pressure peaks are almost twice as high as the negative ones, at the focal point. This indicated that the signal becomes highly nonlinear as it reaches the focal zone.

Figure 4: Acoustic pressure at the water/tissue interface and at the focal point.

Figure 5 shows the frequency content of the recorded signals. The frequency content of the signal entering the tissue phantom is not much different from the source signal. That is, the propagation up to this point is linear. The signal at the focal point is however highly nonlinear, which is seen from the number of generated higher-order harmonics.

The plots of the minimum and maximum on-axis pressure calculated are depicted in Figure 6. The values are computed according to the following formulas

(1)

and are normalized with respect to the maximum pressure at the focal point. Figure 6 shows that p+ is almost twice as high as p-, which is a distinctive feature of nonlinear propagation of a HIFU signal. The bounds of the focal zone are also seen in Figure 6.

Figure 5: Frequency spectrum of the acoustic pressure at the water/tissue interface and at the focal point.

Figure 6: Frequency spectrum of the acoustic pressure at the water/tissue interface and at the focal point.

Notes About the COMSOL Implementation

Mesh and Time Explicit Solver

While solving a wave propagation problem, the used mesh typically has to be fine enough to resolve the frequency content of the signal, that is, the specified number of the higher-order harmonics, N, as discussed in the tutorial Nonlinear Propagation of a Cylindrical Wave — Verification Model. However, a distinctive feature of the model is that the source signal is a pulse and therefore the propagating signal is finite in space. This suggests that a fine mesh is only required in a finite part of the computational domain where the signal is located at a certain time, thus saving a lot of DOFs. This is achievable by enabling the option which can be found in the section of the Time Dependent study step.

The initial mesh built in the model resolves the fundamental harmonic. With the default quartic discretization order, the mesh element size should not exceed 1.5 of the fundamental wavelength. The adaptive mesh refinement ensures that a finer mesh is created to resolve the propagating signal. On the node, select for the Adaptation method and clear the check box. Thus the resulting mesh will always resolve the fundamental harmonic. Then make the capture sharp gradients of the acoustic pressure by typing sqrt(comp1.pr^2+comp1.pz^2).

The Nonlinear Pressure Acoustic, Time Explicit physics is based on dG-FEM and uses an explicit time integration schemes. The solver time step obeys the CFL condition to ensure the stability of the time integration method. That is, Δt ≤ hmin/cmax(p), where hmin is the minimum mesh element size and cmax(p) is the maximum wave propagation speed. In the nonlinear case the speed of sound depends on the acoustic pressure. In this model, the mesh follows the nonlinear propagation of the signal and therefore will have smaller elements around the its sharp fronts and larger elements away from it. It is worth selecting the time-stepping method to and enabling the option available in the section of the time-explicit solver. Figure 7 shows the values of the cell wave time scale elte.wtc variable, that the time step is proportional to. The smaller values are found for the smaller mesh elements in the vicinity of the signal, while the larger values are used in the rest of the computational domain.

Figure 7: Cell wave time scale at times t = 10, 20, 30, and 40 μs.

Minimum and maximum pressure computation

The computation of the minimum and maximum acoustic pressure over time according to Equation 1 is performed using the State Variables functionality available after enabling in the dialog box.

The state variables p_max and p_min are initialized to the acoustic pressure and updated by computing the maximum and minimum between themselves and the acoustic pressure at each time step.

References

1. M.A. Diaz, M.A. Solovchuk, and T.W.H. Sheu, “A conservative numerical scheme for modeling nonlinear acoustic propagation in thermoviscous homegeneous media”, J. Comp. Phys., vol. 363, 2018.

Acoustics_Module/Ultrasound/hifu_tissue_sample

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select .

Click .

Click

In the tree, select .

Click

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file hifu_tissue_sample_parameters.txt.

Rectangle 1 (rect1)

In the toolbar, click

and choose .

In the window for , locate the section.

In the text field, type 0.

In the text field, type 5*T0.

Click to expand the section. Clear the check box.

Create a five-cycle tone burst source signal and inspect its frequency content.

Pulse

In the toolbar, click

and choose .

In the window for , type Pulse in the text field.

In the text field, type pulse.

Locate the section. In the text field, type sin(omega0*t)*(1 - cos(omega0*t/5))*rect1(t).

In the text field, type t.

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

In the text field, type 1.

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


Argument	Lower limit	Upper limit
t	0	9*T0

Click

Results

Tone Burst Pulse

In the window for , type Tone Burst Pulse in the text field.

Function 1

In the window, expand the node, then click .

In the window for , locate the section.

From the list, choose .

Select the check box.

In the text field, type 3*f0.

In the toolbar, click

The plot shows that the frequency content of the signal lies close to the center frequency f0.

Geometry 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Circle 1 (c1)

In the toolbar, click

In the window for , locate the section.

In the text field, type r_source.

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

Click

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type w_source.

In the text field, type r_source-sqrt(r_source^2-w_source^2).

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

Rectangle 2 (r2)

In the toolbar, click

In the window for , locate the section.

In the text field, type r_model+th_abs.

In the text field, type h_model+th_abs-z_tissue.

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	th_abs

Select the check box.

Clear the check box.

Select the check box.

Click

Click the

button in the toolbar.

Rectangle 3 (r3)

In the toolbar, click

In the window for , locate the section.

In the text field, type r_model+th_abs.

In the text field, type z_tissue-(r_source-sqrt(r_source^2-w_source^2)).

Locate the section. In the text field, type r_source-sqrt(r_source^2-w_source^2).

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


Layer name	Thickness (mm)
Layer 1	th_abs

Select the check box.

Clear the check box.

Click

Ignore Edges 1 (ige1)

In the toolbar, click

and choose .

On the object , select Boundary 3 only.

In the toolbar, click

Create two domain probe points on the symmetry axis—one at the water/tissue interface and the other at the focal point—which will record the signal during its propagation from the source.

Definitions

Domain Point Probe 1

In the toolbar, click

and choose .

In the window for , locate the section.

In row , set to z_tissue.

Point Probe Expression 1 (ppb1)

In the window, expand the node, then click .

In the window for , locate the section.

In the text field, type nate.p_t/P0.

Domain Point Probe 2

In the toolbar, click

and choose .

In the window, click .

In the window for , locate the section.

In row , set to F.

Point Probe Expression 2 (ppb2)

In the window, click .

In the window for , locate the section.

In the text field, type nate.p_t/P0.

It is also of interest to compute the minimum and maximum acoustic pressure on the symmetry axis over the simulated time interval.

Click the

button in the toolbar.

In the dialog box, select in the tree.

In the tree, select the check box for the node .

Click .

State Variables 1 (state1)

In the toolbar, click

and choose .

In the window for , locate the section.

In the table, enter the following settings:


State	Initial value	Update expression	Description
p_max	nate.p_t	max(nate.p_t, p_max)	Maximum pressure
p_min	nate.p_t	min(nate.p_t, p_min)	Minimum pressure

Locate the section. From the list, choose .

Select Boundary 2 only.

Now, set up the absorbing layers (sponge layers) used to truncate the computational domain.

Absorbing Layer 1 (ab1)

In the toolbar, click

Select Domains 3–6 only.

In the window for , locate the section.

From the list, choose .

Nonlinear Pressure Acoustics, Time Explicit (nate)

Nonlinear Pressure Acoustics, Time Explicit Model 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Impose the boundary condition upon the interface between the water and the tissue parts of the computational domain.

Material Discontinuity 1

In the toolbar, click

and choose .

Select Boundaries 3 and 11 only.

Pressure 1

In the toolbar, click

and choose .

Select Boundary 18 only.

In the window for , locate the section.

In the p0(t) text field, type P0*pulse(t).

Impedance 1

In the toolbar, click

and choose .

Select Boundaries 6 and 14–17 only.

Set up the materials. The physics interface and the chosen fluid model will suggest which material properties should be defined.

Materials

Water

In the window, under right-click and choose .

In the window for , type Water in the text field.

Select Domains 1 and 4 only.

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


Property	Variable	Value	Unit	Property group
Density	rho	rho_water	kg/m³	Basic
Speed of sound	c	c_water	m/s	Basic
Sound diffusivity	delta_diff	delta_water	m²/s	Attenuation and dissipation model
Parameter of nonlinearity	BA	BA_water	1	Nonlinear model

Tissue

Right-click and choose .

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

Select Domains 2, 3, 5, and 6 only.

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


Property	Variable	Value	Unit	Property group
Density	rho	rho_tissue	kg/m³	Basic
Speed of sound	c	c_tissue	m/s	Basic
Sound diffusivity	delta_diff	delta_tissue	m²/s	Attenuation and dissipation model
Parameter of nonlinearity	BA	BA_tissue	1	Nonlinear model

Now, proceed to the mesh. Even though the problem at hand is nonlinear, make the mesh resolve only the fundamental harmonic. To do it, limit the maximum element size by 1.5 of the wavelength in the corresponding materials.

Mesh 1

Free Triangular 1

In the toolbar, click

Size 1

Right-click and choose .

In the window for , locate the section.

From the list, choose .

Select Domains 1 and 4 only.

Locate the section. Click the button.

Locate the section. Select the check box.

In the associated text field, type c_water/f0/1.5.

Size 2

In the window, right-click and choose .

In the window for , locate the section.

From the list, choose .

Select Domains 2, 3, 5, and 6 only.

Locate the section. Click the button.

Locate the section. Select the check box.

In the associated text field, type c_tissue/f0/1.5.

Click

Study 1

Step 1: Time Dependent

In the window, under click .

In the window for , locate the section.

In the text field, type range(0, T0, 55*T0).

This setting saves the solution at times multiple to T0 in the whole computational domain. It only influences the stored solution (and thus the file size). The internal time steps taken by the solver are automatically controlled by COMSOL to fulfill the appropriate CFL condition.

On the other hand, the signals at the probes will be computed for each time step taken by the solver thus providing a much higher temporal resolution of the results.

Click to expand the section. Select the check box.

This setting enables the adaptation of the mesh with respect to the propagating signal by creating intermediate meshes. Those will be fine enough in the parts of the computational domain where it is required to resolve higher-order harmonics generated because of the nonlinear nature of the problem.

Solution 1 (sol1)

In the toolbar, click