COMSOL 6.4 - RF Implant Heating in MRI

Magnetic resonance imaging (MRI) relies on the application of time-varying magnetic fields. Among those, radio frequency (RF) fields are used for volumetric excitation of the target. This tutorial model focuses on the RF excitation and its heating effect on a passive conductive implant inserted into anatomical vertebrae¹ and plunged into a jelly phantom having human body’s tissue properties.

RF-induced heating occurs any time a material with finite mass is invested by RF power; this is usually referred to as specific absorption rate (SAR):

(1)

where σ is the conductivity of material, ρ is the mass density, and |E| is the norm of the electric field (RMS). This means that a human body’s tissue yields a finite SAR even without a conductive prosthesis. However, when a passive conductive implant is present in the tissue, a secondary RF field is induced in the implant that in turn generates an additional, dominant (compared to the original excitation), scattered field toward the surroundings. As a result, the RF power absorption by the tissue, and thus the SAR, are locally enhanced. The resulting local heating increase and consequent tissue damage represent one of the major hazards for the patient under MRI test.

This MRI implant heating tutorial model shows the RF-induced heating on a passive conductive polyaxial screw implant inserted into C3-C4-C5 anatomical vertebrae and plunged into a jelly phantom brick so as to mimic the tissue enclosing the spine. The phantom has standard physical properties and geometry as reported in the ASTM‑F2182‑19e2 standard. The SAR is computed in a first Frequency Domain study step. The resulting heating is coupled to a second Time Dependent study step to assess temperature increase consistently with the procedure reported in the ASTM‑F2182-19e2 standard.

Model Definition

The geometry of the system is shown in Figure 1. The 3D model consists of a birdcage coil for RF excitation, a phantom brick having human body’s tissue properties, a series of C3-C4-C5 vertebrae, and a polyaxial screw implant. The vertebreae and polyaxial screws are taken directly from the tutorial STL Import 2 — Combining Geometry with an Imported Mesh in the COMSOL Multiphysics Application Library; they are shown enlarged in Figure 2 for better visualization.

The RF excitation is the same used in MRI coil tutorial model; shortly: two lumped ports are used for quadrature excitation and a number of lumped-element capacitors are included for fine frequency tuning and field uniformity inside the MRI chamber.

The 3D mesh is shown in Figure 3. It is built using a free tetrahedral mesh for all domains except for a swept mesh with a distributed number of elements in the outer PML domain.

The implementation requires the use of the Electromagnetic Waves, Frequency Domain interface and the Heat Transfer in Solids interface. The former interface is used to compute the SAR in the phantom and vertebrae at 63.87 MHz operation frequency. The resulting heating is then incorporated into the latter interface via a multiphysics coupling in order to evaluate the temperature evolution over 6 minutes of steady RF-heating exposure.

Figure 1: 3D geometry of the MRI system including a birdcage coil, a phantom brick, three anatomical vertebrae and a polyaxial screw implant.

Figure 2: Zoomed view on the polyaxial screw implant inserted into C3-C4-C5 vertebrae.

Figure 3: Computational mesh used in this tutorial.

Results and Discussion

Two study steps are performed in this tutorial: a Frequency Domain study step to evaluate the SAR and a Time Dependent study step to assess the temperature rise consequent to the RF heating.

Figure 4 displays, on the primary y-axis, the temperature rise over time (with respect to the baseline) at two point probes located at the distal ends of the implant’s rod. The secondary y-axis reports the background SAR at a reference point probe in the phantom specular to the implant location.

Figure 5 shows the electric field norm on the implant surface remarking that the field strengthen at the distal ends of the implant. This field enhancement stems from the scattering induced by the implant.

Figure 6 depicts the SAR on a phantom region surrounding the vertebrae and the implant. The SAR is maximum at the same locations where the RF field strengthen (see Figure 5).

Figure 7 shows the temperature distribution on the implant after 6 minutes of steady RF-heating exposure. This result is consistent with the plot of Figure 4.

Figure 4: (left axis) Temperature evolution at rod probes, (right axis) SAR at reference probe.

Figure 5: Electric field norm on the implant surface.

Figure 6: SAR on a phantom region surrounding the implant and the vertebrae.

Figure 7: Temperature distribution on the implant after 6 min of RF-heating exposure.

Reference

1. Reprinted, with permission from ASTM‑F2182-19e2, Standard Test Method for Measurement of Radio Frequency Induced Heating On or Near Passive Implants During Magnetic Resonance Imaging, copyright ASTM International. A copy of the complete standard may be obtained from ASTM, www.astm.org.

RF_Module/Microwave_Heating/mri_implant_heating

Modeling Instructions

From the menu, choose .

In the window, click

In the window, click

From the menu, choose .

Browse to the model’s Application Libraries folder and double-click the file stl_2_combine_geom_mesh_geom_sequence.mph.

Global Definitions

Parameters - Vertebrae and Screw

In the window, under click .

In the window for , type Parameters - Vertebrae and Screw in the text field.

Parameters - Birdcage Coil

In the toolbar, click

and choose > .

In the window for , type Parameters - Birdcage Coil in the text field.

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


Name	Expression	Value	Description
c_value	29.5[pF]	2.95E-11 F	Capacitance used on the rungs
r_coil	xPhantom		Radius of the coil
h_coil	1.5*zPhantom		Height of the coil
l_element	0.01[m]	0.01 m	Length of the capacitive elements
t_ring	0.015[m]	0.015 m	Thickness of the ring
f0	63.87[MHz]	6.387E7 Hz	Operating frequency
wl0	c_const/f0	4.6938 m	Free space wavelength
rAir	1.2[m]	1.2 m	Radius of air sphere
tPml	0.2[m]	0.2 m	PML thickness
V0	540[V]	540 V	Voltage excitation

Parameters - Phantom

In the toolbar, click

and choose > .

In the window for , type Parameters - Phantom in the text field.

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


Name	Expression	Value	Description
xPhantom	42[cm]	0.42 m	Phantom length along x-axis
yPhantom	9[cm]	0.09 m	Phantom length along y-axis
zPhantom	65[cm]	0.65 m	Phantom length along z-axis
epsrPhantom	80	80	Relative permittivity phantom
sigmaPhantom	0.47[S/m]	0.47 S/m	Electrical conductivity phantom
murPhantom	1	1	Relative permeability phantom
cPhantom	4150[J/kg/K]	4150 J/(kg·K)	Heat capacity phantom
alphaPhantom	1.3e-7[m^2/s]	1.3E-7 m²/s	Thermal diffusivity phantom
volPhantom	xPhantomyPhantomzPhantom	0.02457 m³	Volume of phantom
rhoNaCl	1.32[g/L]	1.32 kg/m³	Density of NaCl
rhoPAA	10[g/L]	10 kg/m³	Density of polyacrylic acid (PAA)
rhoPhantom	rhoWater+ rhoNaCl + rhoPAA	1004.6 kg/m³	Density of phantom
rhoWater	993.31[kg/m^3]	993.31 kg/m³	Density of water at 37C
kPhantom	alphaPhantomrhoPhantomcPhantom	0.542 W/(m·K)	Thermal conductivity of phantom
dWall	2[cm]	0.02 m	Distance from phantom wall

First, build the geometry of the polyaxial screw and rod.

Geometry 1

Delete Entities 1 (del1)

In the window, expand the > node, then click .

In the window for , click

Second, import the geometry sequence for the birdcage coil and phantom, then build all.

In the toolbar, click and choose .

Browse to the model’s Application Libraries folder and double-click the file mri_implant_heating_geom_sequence.mph.

In the toolbar, click

Now define useful selections under the Geometry node.

Non-PML

In the toolbar, click

and choose .

In the window for , type Non-PML in the text field.

Locate the section. In the text field, type -0.705[cm].

In the text field, type -xPhantom/2+dWall+d_rod/2+1.15*c3_zw.

In the text field, type -(c3_xw+c5_xw)/2.

Locate the section. In the text field, type rAir-tPml/2.

Locate the section. From the list, choose .

Click

PML

In the toolbar, click

and choose .

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

Locate the section. Click

In the dialog, select in the list.

Click .

In the window for , click

Non-PML Boundaries

In the window, under > right-click and choose .

In the window for , type Non-PML Boundaries in the text field.

Locate the section. From the list, choose .

Build the mesh-based geometry and perform a mesh cleanup.

Mesh-Based Geometry 1

Union 1

In the window, expand the > node, then click .

In the window for , locate the section.

From the list, choose .

In the text field, type 0.01[mm].

In the window, right-click and choose .

Right-click and choose > .

Cleanup Wizard

Go to the window.

Click the button.

In the text field, type 0.1[mm].

Click the button in the window toolbar.

Definitions (comp1)

Step 1 (step1)

In the toolbar, click

and choose .

In the window for , locate the section.

In the text field, type 0.25.

Continue setting selections under the Definitions node.

Vertebrae and Screws

In the toolbar, click

In the window for , type Vertebrae and Screws in the text field.

Locate the section. Under , click

In the dialog, in the list, choose and .

Click .

Discs

In the toolbar, click

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

Locate the section. Under , click

In the dialog, select in the list.

Click .

Vertebrae

In the toolbar, click

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

Locate the section. Under , click

In the dialog, select in the list.

Click .

In the window for , locate the section.

Under , click

In the dialog, select in the list.

Click .

Port 1

In the toolbar, click

In the window for , type Port 1 in the text field.

Locate the section. From the list, choose .

Click

In the dialog, type 36 in the text field.

Click .

Port 2

In the toolbar, click

In the window for , type Port 2 in the text field.

Locate the section. From the list, choose .

Click

In the dialog, type 21 in the text field.

Click .

Phantom

In the toolbar, click

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

Select Domain 6 only.

HT Domain

In the toolbar, click

In the window for , type HT Domain in the text field.

Locate the section. Under , click

In the dialog, in the list, choose , , , and .

Click .

Implant

In the toolbar, click

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

Locate the section. Under , click

In the dialog, in the list, choose and .

Click .

Implant Boundary

In the toolbar, click

In the window for , type Implant Boundary in the text field.

Locate the section. Under , click

In the dialog, select in the list.

Click .

SAR Domain

In the toolbar, click

In the window for , type SAR Domain in the text field.

Locate the section. Under , click

In the dialog, select in the list.

Click .

In the window for , locate the section.

Under , click

In the dialog, select in the list.

Click .

All Elements

In the toolbar, click

In the window for , type All Elements in the text field.

Locate the section. From the list, choose .

Click

In the dialog, type 19, 21, 23, 24, 26, 28, 29, 31, 33, 34, 36, 38, 171, 173, 175, 176, 178, 180, 181, 183, 185, 186, 188, 190 in the text field.

Click .

Temperature Probes

In the toolbar, click

In the window for , type Temperature Probes in the text field.

Locate the section. From the list, choose .

Select Points 71 and 84 only.

SAR Probe

In the toolbar, click

In the window for , type SAR Probe in the text field.

Locate the section. From the list, choose .

Select Point 102 only.

Now add the following physics interfaces with the corresponding settings.

Add Physics

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 tree, select > .

Click the button in the window toolbar.

In the text field, type event.

In the tree, select > > .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Electromagnetic Waves, Frequency Domain (emw)

Perfect Electric Conductor 2

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Perfect Electric Conductor 3

In the toolbar, click

and choose .

In the window for , locate the section.

Click

In the dialog, type 13, 14, 103, 116 in the text field.

Click .

Lumped Port 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

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

Lumped Port 2

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

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

In the θin text field, type pi/2.

Lumped Element 1

In the toolbar, click

and choose .

Select Boundaries 19, 23, 24, 26, 28, 29, 31, 33, 34, 38, 171, 173, 175, 176, 178, 180, 181, 183, 185, 186, 188, and 190 only.

In the window for , locate the section.

From the list, choose .

In the Celement text field, type c_value.

Click the button in the window toolbar.

Specific Absorption Rate 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Impedance Boundary Condition 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Definitions (comp1)

Perfectly Matched Layer 1 (pml1)

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

Heat Transfer in Solids (ht)

In the window, under click .

In the window for , locate the section.

From the list, choose .

Heat Flux 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

In the h text field, type 5.

Locate the section. From the list, choose .

Heat Source 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the Q0 text field, type emw.Qrh*switch*step1(t/1[s]).

Boundary Heat Source 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Locate the section. In the Qb text field, type emw.Qsrh*switch*step1(t/1[s]).

Events (ev)

In the window, under click .

Discrete States 1

In the toolbar, click

and choose .

In the window for , locate the section.

In the table, enter the following settings:


Name	Initial value (u0)	Description
switch	1

Explicit Event 1

In the toolbar, click

and choose .

In the window for , locate the section.

In the ti text field, type 360[s].

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


	Expression
switch	0

At this point, add and assign the following materials to the respective domains. Set relative permeability to one for all materials.

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.

Materials

Phantom Material

In the window, under right-click and choose .

In the window for , type Phantom Material in the text field.

Select Domain 6 only.

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


Property	Variable	Value	Unit	Property group
Relative permittivity	epsilonr_iso ; epsilonrii = epsilonr_iso, epsilonrij = 0	epsrPhantom	1	Basic
Relative permeability	mur_iso ; murii = mur_iso, murij = 0	murPhantom	1	Basic
Electric conductivity	sigma_iso ; sigmaii = sigma_iso, sigmaij = 0	sigmaPhantom	S/m	Basic
Density	rho	rhoPhantom	kg/m³	Basic
Thermal conductivity	k_iso ; kii = k_iso, kij = 0	kPhantom	W/(m·K)	Basic
Heat capacity at constant pressure	Cp	cPhantom	J/(kg·K)	Basic

Add Material

Go to the window.

In the tree, select > > > > .

Click the button in the window toolbar.

Materials

Ti - 6 Al - 4 V (Grade 5) [solid,annealed] (mat3)

In the window for , locate the section.

From the list, choose .

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


Property	Variable	Value	Unit	Property group
Relative permittivity	epsilonr_iso ; epsilonrii = epsilonr_iso, epsilonrij = 0	1	1	Basic
Relative permeability	mur_iso ; murii = mur_iso, murij = 0	1	1	Basic

Add Material

Go to the window.

In the tree, select > > .

Click the button in the window toolbar.

Materials

Bone Cancellous (mat4)

In the window for , locate the section.

From the list, choose .

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


Property	Variable	Value	Unit	Property group
Relative permeability	mur_iso ; murii = mur_iso, murij = 0	1	1	Basic

Add Material

Go to the window.

In the tree, select > .

Click the button in the window toolbar.

Materials

Bone (mat5)

In the window, under > right-click and choose .

Complete missing parameters for the Cartilage material using those of the phantom and assign it to the discs between vertebrae.

Add Material

Go to the window.

In the tree, select > > .

Click the button in the window toolbar.

In the toolbar, click

to close the window.

Materials

Cartilage (mat5)

In the window for , locate the section.

From the list, choose .

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


Property	Variable	Value	Unit	Property group
Relative permeability	mur_iso ; murii = mur_iso, murij = 0	1	1	Basic
Density	rho	rhoPhantom	kg/m³	Basic
Thermal conductivity	k_iso ; kii = k_iso, kij = 0	kPhantom	W/(m·K)	Basic
Heat capacity at constant pressure	Cp	cPhantom	J/(kg·K)	Basic