Development of a Magnetic Fluid Heating FEM Simulation Model with Coupled Steady State Magnetic and Transient Thermal Calculation

: Magnetic ﬂuid hyperthermia has gained much attention in recent years due to its potential in cancer treatment. Magnetic ﬂuid is a colloidal liquid made of nanoscale magnetic particles suspended in a carrier ﬂuid. The properties of a commercial magnetic ﬂuid consisting of maghemite ( γ -Fe 2 O 3 ) particles suspended in mineral oil were used in the scope of our research. The paper deals with a novel approach to the development of a magnetic ﬂuid FEM model of a laboratory setup, with consideration of the electromagnetic steady state and thermal transient calculation soft coupling. Also, adjustment of the mathematical model was added in such a way that it enables a link between the magnetic and thermal calculations in commercial software. The effective anisotropy’s inﬂuence on the calculations is considered. The simulation was done for different magnetic ﬁeld parameters. The initial temperature was also varied so that a direct comparison could be made between the simulation and the measurements. A good indicator of the accuracy of the simulation are the SAR values. The relative differences in SAR values were in the range from 4.2–24.9%. Such a model can be used for assessing the heating performance of a magnetic ﬂuid with selected parameters. It can also be used to search for the optimal parameters required to design an optimal magnetic ﬂuid.


Introduction
Magnetic fluid hyperthermia has gained much attention in recent years due to its potential in cancer treatment. Hyperthermia treatment should impair tumour cells without damaging the healthy surrounding tissue. In practice, the treatment involves heating the tumour region to the temperature range of 40-44 • C [1]. Therefore, it is essential to know the heating characteristics of the magnetic fluid. In addition to heating of the cells, it should be mentioned that there are other possible methods, like freezing the cells [2].
Magnetic fluid, also named ferrofluid, is a colloidal liquid made of nanoscale magnetic particles suspended in a carrier fluid. The particles are coated with a surfactant to prevent agglomeration. The properties of a commercial magnetic fluid consisting of maghemite (γ-Fe 2 O 3 ) particles suspended in mineral oil were used in the scope of our research.
Experimentally, magnetic fluid heating characteristics are often measured using a method called the calorimetric method. There, the magnetic fluid is subjected to an alternating magnetic field, and the time dependence of the temperature is measured, as has been done in [3]. The heating characteristics of the magnetic fluids are usually represented with a parameter called the Specific Absorption Rate (SAR) in Watts per gram of magnetic nanoparticles (MNPs). The value of the SAR can be obtained directly from the time-temperature curve. The most straightforward approach in determining SAR is using the initial slope method, where the SAR is calculated from the initial slope of the time-temperature curve as: Here, c f is the fluid's specific heat, ρ f is the density of the fluid, m MNP is the mass of the MNPs in the fluid and (∆T/∆t) i is the initial slope of the curve. This approach assumes that the heating process is adiabatic, or that cooling is negligible. It was used in the literature [1,[4][5][6].
Due to the difficulties and time consumption in characterizing magnetic fluids experimentally, a numerical model was developed using the Finite Element Method (FEM). With this, it is possible to simulate the heating of the magnetic fluid in a calorimetric experiment environment and obtain results for magnetic fluids with different parameters much faster. It can also be used for the determination of magnetic fluid parameters for optimal heating.
Related work can be found in the literature. In the literature [7], the authors deal with the development of the analytical relationship and computations of power dissipation in magnetic fluid in an alternating magnetic field based on relaxation processes in the magnetic fluid. Similarly, article [8] presents a model based on relaxation processes to calculate magnetic losses and SAR based on the parameters of the magnetic fluid. The model was validated with the use of calorimetric measurements for different magnetic field amplitudes, frequencies and shapes. Article [9] deals with the determination of the temperature dependence of the heating power of the magnetic fluid. Temperature dependence was determined experimentally via magnetic measurements. The obtained curves were used in transient FEM analysis (using OPERA software) of the heating power. The results were compared with the measurements. Literature [10] deals with the development of a numerical model for estimating temperature distribution in a magnetic fluid system, validating it by analysis of the temperature distribution between experimental measurements and numerical analysis by FEM (using ABAQUS software, ABAQUS Inc., Palo Alto, CA, USA). Initially, experiments were conducted and the temperature increase was measured. For the obtained curves, Specific Loss Power (SLP) was determined and used for numerical modelling of the system. Many articles deal with numerical simulations calculating the heating of MNPs in cancer tissue using a heating model based on relaxation processes and Pennes bioheat equation [11][12][13][14][15][16]. Similarly, [17] deals with the numerical modelling of hypethermia using the Pennes model for heat transfer while using the Gaussian heat generation distribution proposed by [18]. In [19], heating of the MNPs was simulated using a stochastic Néel-Brown Langevin equation and Monte Carlo simulations. The results were compared with the calorimetric measurements results. In [20], a 2D axisymmetric model was developed in COMSOL using heat transfer equations and Navier-Stokes equations. For heating power, a temperature-dependent function was used, obtained via an optimization study. In [21], the efficiency of MNPs' heating was investigated in the dependence of various magnetic fluid parameters, including particle size, shape, anisotropy, and the degree of aggregation or agglomeration of the particles. The aggregation of particles in magnetic fluids tends to be avoided, but can be useful in some treatment applications, like MRI-guided drug delivery, as demonstrated in [22][23][24]. While simulations of the heating of MNPs inside a tumour give important results, they do not, however, describe the intrinsic heating characteristic of the magnetic fluid. In the literature, mathematical models used to calculate the heating of the magnetic fluid are known. They are used to simulate heating of the MNPs in cancer tissue in [11][12][13][14][15][16]. For the laboratory setting, simulations in the literature use pre-measured heating results for heating power data, and simulate temperature distribution using the transient thermal calculations in [9,10], although to our knowledge, coupled magnetic and thermal simulations for the laboratory setting are missing. Therefore, our work focuses on the heating capabilities of a magnetic fluid in a laboratory setting using coupled steady state magnetic and transient thermal calculations for the simulation.
The main contributions of our work are: • Development of a magnetic fluid FEM model with consideration of electromagnetic and thermal soft coupling.

•
Adjustment of the mathematical model in such a way that it enables a link between the magnetic and thermal calculations in commercial software (Altair Flux3D).

•
Analysis of an effective anisotropy's influence on the calculation results, and suggestions for correct setting of an effective anisotropy.

•
The correctness of the model and the correctness of the used parameters are confirmed with measurements.
The paper consists of seven Sections. The measurement system for measuring the heating characteristics of the magnetic fluid is described in Section 2. The mathematical model used for calculations is described in Section 3. The simulation model and simulation process are described in Section 4. The parameters for the mathematical model and material properties used in the simulation are described in Section 5 in Sections 5.1 and 5.2, respectively. The results are gathered in Section 6, and conclusions are in Section 7.

Measurement System
The measurement system in our laboratory was used for the measuring (Laboratory for Applied Electromagnetics, Faculty of Electrical Engineering and Computer Science, University of Maribor). The measuring system is based on the calorimetric method for measuring magnetic fluids. It consists of a supply coil, a glass container with magnetic fluid, two field measuring coils and an optical temperature sensor. The supply coil is cooled with water. A function generator was used to generate the appropriate magnetic field, whose signal was then amplified with an amplifier. Additional capacitors were used to create an LC resonator with the supply coil. The setup is presented schematically in Figure 1. The measurement procedure is as follows. First, a sinusoid signal is generated with the function generator. The signal is then amplified using a power amplifier. Beforehand, a suitable combination of capacitors is selected, to get the capacitance C that, combined with the supply coil's inductance L, forms an LC resonator with the desired frequency f. This is the frequency set on the function generator. The supply coil generates an alternating magnetic field that heats the magnetic fluid positioned within the glass container inside the supply coil, as shown in Figure 1. The temperature is measured with an optical thermometer and the data are saved to the computer. Two measurement coils are used to measure the magnetic field in the magnetic fluid. The process of field measurement is described in [3] in more detail. From the obtained time-temperature curve at the measured magnetic field, SAR is calculated according to Equation (1).

Mathematical Model
To calculate the heating power of the magnetic fluid, a model based on the relaxation of magnetization was used, as described in [7,8]. There are two mechanisms of relaxation. One is relaxation by physical rotation of the MNP, known as Brownian relaxation. It is described by the relaxation time: where V h is the hydrodynamic volume of the MNP, which takes into account the surfactant layer on top of the MNP, η is the viscosity of the carrier fluid, k B is the Boltzmann constant and T is the temperature in kelvin. The other type of relaxation is the flipping of the magnetization of the MNP, known as Néelian relaxation. It is described by the relaxation time: where τ 0 is the attempt time, K eff is the effective anisotropy constant of the MNP and V p is the volume of the MNP. The effective relaxation time τ is given as: Volumetric power dissipation in units of W/m 3 is expressed as: as was shown by Rosensweig in [7]. Here, µ 0 is the permeability of free space, f frequency, H magnetic field strength amplitude and χ the imaginary part of the complex susceptibility, expressed as: where ω = 2πf is the angular frequency and χ 0 is the equilibrium susceptibility: with χ i being initial susceptibility and L(α) the Langevin function and α the Langevin parameter. They are expressed as: respectively. M s represents the saturation magnetization of the bulk maghemite, and φ represents the volume fraction of MNPs.
To take into account the polydisperse nature of MNP sizes, the particles were divided into size groups. The heating power is the sum of the contributions of individual groups: Here, index j represents the individual group and M is the number of groups. The size-dependent variables contained in χ j are: Particle volume V pj and hydrodynamic volume V hj , • Brown and Néel relaxation times τ Bj and τ Nj , • Langevin parameter α j .

Simulation
The Altair Flux 2018.1.2 software (Altair Engineering Inc., Troy, MI, USA) was used for the FEM simulation. The Flux3D was used for modelling, and Steady State AC Magnetic coupled with Transient Thermal option was used for solving, with the adjusted mathematical model used to couple them. The model that was built is presented schematically in Figure 2. For the simulation, based on previous information, a glass container was modelled with a marked height of the magnetic fluid.
A function offered by the Flux software named Infinite box was used for the boundary conditions. The technique uses a transformation to model an infinite domain. In our model, the infinite box was made by using two superimposed parallelepipeds. The faces of the outer parallelepiped represent the images of the infinity, where the magnetic field and its potential are equal to zero. The temperature at infinity is equal to the initial temperature of the problem. Afterwards, faces and volumes were built, and the problem was meshed. A circular non-meshed coil of a rectangular cross-section was used for the supply coil, with the dimensions shown in Figure 2.
The regions were defined as the following: A non-active, or air region, in the infinity box region and the region between the infinity box and the container, glass for the container, magnetic fluid in the container below the magnetic fluid height marking, and air in the container above the marking. It is important to define the appropriate area for radiation and convection heat exchange to get correct results. For the presented problem, a face region with radiation and convection heat exchange was defined on the outer boundary of the container and on the upper face of the air inside the container. The regions are represented in Figure 3. Equations (2)-(11) of the mathematical model were defined as spatial quantities in the scope of the Flux3D software environment. The model works as follows. The programme calculates the steady state magnetic conditions, and, after that, the heating power in the magnetic fluid, using Equations (2)-(11), initial temperature, magnetic field strength amplitude and frequency in the first time step. Based on the power, it calculates the temperature for the next time step. Then the process repeats for the new time step until the end of the simulation. The simulation flowchart is presented in Figure 4. After the simulation ended, the average magnetic field strength amplitude in the whole magnetic fluid volume was calculated for each time step and the average value over time was taken. To get the temperature and SAR values, one point inside the magnetic fluid (the position of a temperature measuring sensor used in the case of measurement) was chosen to simulate the temperature measurement with a temperature probe. The location of the reading is shown in Figure 5, together with the Infinite box boundary conditions. To get the value of SAR, we used equation [8]: Here, ρ MNP represents the mass density of MNPs in the magnetic fluid, and is equal to ρ MNP = φρ m , where ρ m = 4.9 g/cm 3 is the density of the bulk maghemite. Equation (12) gives us a temperature-dependent function of SAR, the value of which at t = 0 s is directly comparable to the result of Equation (1).
The results given by this model were validated using the experimental measurement results obtained by the measurement system described in Section 2. It was done by comparing the time-temperature curves and initial SAR values from the simulations and measurements.

Selection of Parameters
In the modelling of the simulation, careful considerations must be made regarding parameters. Two types of parameters were present-the ones used in the mathematical model for the heating power calculation and the ones describing the magnetic and thermal properties of the materials.

Parameters Used for the Determination of the Software's Spatial Quantities Determined with the Equations
In the mathematical model, the essential parameters are the volume fraction of MNPs φ with a corresponding size distribution, volumes V p and V h , the viscosity of carrier fluid η, attempt time τ 0 , effective anisotropy K eff and the magnetization of the bulk material M s .

Volume Fraction φ
Volume fraction φ represents the content of MNPs in the magnetic fluid. For our sample of magnetic fluid, the volume fraction was measured at φ = 10.57%. Because the particles in the magnetic fluid are polydispersed, the particles were divided into ten groups of different sizes. The sizes and volume fractions are presented in Table 1 and Figure 6.

Volumes V p and V h
Volume V p is the volume of one MNP. For the volume calculation, the particles were assumed to be spherical. Therefore, V p can be estimated as V p = π/6d 3 , where d is the diameter of the MNP. Volume V h includes the surfactant layer s on the particle to the volume. The surfactant layer was estimated to be s = 1 nm, which adds 2 nm to the diameter of the particle.

Viscosity η
The following function was fitted to the data of temperature dependence of viscosity for mineral oil [26]: Here, the fitting constants are: A = 3.4313 × 10 −5 Pa s, B = 1248 K, C = −155.05 K.

Attempt Time τ 0
The attempt time generally depends on the MNP material properties and, according to the literature [27], falls between 10 −13 and 10 −8 s. The value of 10 −9 s was assumed, the same as in [8].

Effective Anisotropy K eff
This is one of the most critical parameters and one of the hardest to determine. Magnetic anisotropy describes the dependence of the magnetic properties of the material on the direction. It can be defined as the energy required to rotate the magnetization from the easy to the hard axis [28]. Anisotropy results from several contributions: Magnetocrystalline, shape and surface are some of them. The magnetic moment of the MNP generally coincides with the axis of the magnetic anisotropy [29]. Together with the MNP volume, the effective anisotropy coefficient describes the energy barrier the magnetization has to overcome to flip the direction (energy E = K eff V).
In the literature the measurements of maghemite MNPs effective anisotropy yielded different results in the range from 10 3 to 10 4 J/m 3 [30][31][32]. To determine the value of K eff , the Stoner-Wolfarth model for shape anisotropy was used, described by equation (14) [33]. Using a Transmission Electron Microscope (TEM) picture, the aspect ratios r of the MNPs was determined to fall in the range between 1 and 2 for our sample, with the average being around 1.4.
Here, N t and N l are transverse and longitudinal demagnetizing factors, respectively. They are related by N l + 2N t = 2π. The longitudinal demagnetizing factor of a spheroid is calculated as: where = √ (1 − 1/r 2 ) is the eccentricity. Combining shape anisotropy with magnetocrystalline anisotropy, which is 4.6 kJ/m 3 for maghemite [7], the upper limit for the value of effective anisotropy becomes 17.3 kJ/m 3 , which was the value used in the simulation.

Saturation Magnetization of the Bulk Material M s
As the name suggests, the parameter M s represents the saturation magnetization of the bulk magnetic material. In our case, that was maghemite. The saturation magnetization of bulk maghemite is 400 kA/m [8].
The important parameters used for spatial quantities are collected in Table 2. Table 2. Values of parameters (for spatial quantities) used in the simulation.

Parameter Value
Attempt time τ 0 10 −9 s Effective anisotropy K eff 17.3 kJ/m 3 Saturation magnetization of the bulk material M s 400 kA/m

Basic Parameters Used in the Simulation Model for Material's Properties' Determination
The following material parameters were considered in the simulation model: Density, specific heat capacity, thermal conductivity or radiation and convection coefficients, and the magnetic permeability of glass for the container, air for the surrounding area, and magnetic fluid for the sample in the container. The parameters are: For air:   • For the magnetic permeability, the initial susceptibility was calculated using Equation (8) and converted into permeability using the equation: • Initial magnetic permeability was then ( Figure 9): • Saturation magnetization of the magnetic fluid: M f = φM s = 42.3 kA/m, • The thermal conductivity was calculated with equation [38]: where λ f is the thermal conductivity of the magnetic fluid, and β = (λ p − λ f )/(λ p + 2λ f ), with λ p being the thermal conductivity of the MNP. For the thermal conductivity of the carrier fluid, the value of λ f = 0.15 W/mK was estimated for mineral oil, taking literature [39] as a guide. For maghemite, the thermal conductivity falls in the range of 0.86-1.30 W/mK, according to [40]. For the value λ p = 0.92 W/mK, using Equation (17), the magnetic fluid thermal conductivity was λ = 0.187 W/mK. On the face region ( Figure 3) around the container, the emissivity (value of 0.93) and absorption (value of 0.2) coefficients (for radiation), and convection coefficient (value of 10 W/m 2 K) were defined. These were used to simulate radiation and natural convection losses.

Results
The simulation was run for different magnetic field parameters (field strength and frequency). Also, the initial temperature was varied, so that a direct comparison could be done between the simulation results and the measurements.

Mesh Dependency Test
To test the mesh, a simulation was run for three different meshes at magnetic field parameters of H = 7.83 kHz and f = 281 kHz. Figure 12 shows the comparison of the results for these three meshes. The original mesh was the mesh used for all the other calculations, Rarer mesh has less elements than Original mesh, and Denser mesh has more elements. The differences in the results with different mesh densities used were not significant (less than 0.5 %). The number of elements for used mesh densities are shown in Table 3. Table 3. Number of nodes, line elements, surface elements and volume elements for the three meshes used in the comparison.  Figures 13-15 represent the time evolution of temperature at different conditions (frequency, magnetic field amplitude). Since the sample volume was not the same in measurement and simulation due to the difficult exact determination of the magnetic fluid volume used for the measurement, a deviation in temperature between the two was expected. Although this was the case, the initial slope of both should be comparable and, with it, comparable SAR.   The temperature dependency on time, presented in Figures 13-15, shows a similar dependency as shown in other works, like in the literature [10,16].

Comparison with the Measurements
The SAR values obtained with Equation (1) for measurement and the initial SAR values obtained with Equation (12) for simulation were compared. The results are shown in Table 4. The SAR values are a good indicator of the accuracy of the simulation. Relative differences in SAR values are in the range from 4.2-24.9%.

Effective Anisotropy's Influence on the Simulation Results
A comparison was made for the results of calculations for different values of K eff . The comparison of time-temperature curves and time-SAR curves is shown in Figures 16 and 17, respectively.  In Figure 16, it can be seen that the temperature curves rise with the rising value of K eff to it being around 25 kJ/m 3 and 30 kJ/m 3 (curves (c) and (d) in Figure 16), for the magnetic field with parameters H = 7.83 kA/m and f = 281 kHz. At higher values of K eff , the temperature curves start falling. Looking at the time-SAR curves in Figure 17, it can be seen that time dependence changes with the changing of the K eff . For the given magnetic field, the SAR curves represented in Figure 17  The reason for this behaviour comes down to the dependence of heating power P (SAR and P are linearly dependent) on temperature and K eff , based on Equation (11). This dependence is presented graphically in Figure 18 for the magnetic field parameters of f = 281 kHz and H = 7.83 kA/m. For lower K eff , the power output fell with the temperature increase. For higher values of K eff (around 45 kJ/m 3 and up), the maximum of the power output appeared at a temperature of about 215 • C. It can be seen that the results are heaviliy dependent on the effective anisotropy constant. In the literature [21] ( Figure 5 in the literature [21]), it was shown that heating power reaches a maximum at different particle sizes for different anisotropy constants. This explains the local maximums seen in Figure 18 (heating power dependency on K eff ) at low temperatures.

Magnetic Field Strength Amplitude's Influence on the Simulation Results
Time-temperature curves at different magnetic field strength amplitudes are represented on Figure 19. It can be seen that, with the rising magnetic field strength amplitude, the curves rise, meaning that the heating is more intense at higher magnetic field strengths. This agrees with the findings in the literature [8,20]. Figure 20 shows the initial SAR dependence on magnetic field strength amplitude. As seen in Figure 20, the calculated initial SAR is rising with the magnetic field strength amplitude. This agrees with the results in [8], where the results for alternating magnetic field show a similar SAR-H dependence.

Frequency's Influence on the Simulation Results
According to [8], the heating power should be increasing with the frequency. This is confirmed with Figure 21 which shows time-temperature curves calculated at different frequencies, and Figure 22 which shows how initial SAR rises with frequency.

Discussion of the Results' Accuracy
The accuracy of the results is limited, as is seen from the simulation results deviations from the measurements (4.2-24.9%). Comparing these results with the results from studies that use measured heating power as a simulation input [9,10], the second ones are certainly more accurate. However, they are limited to one specific fluid, while our method could be expanded to magnetic fluids with different parameters. The studies simulating the heating of a tissue with embedded MNPs [11][12][13][14][15][16] give better understanding of the heating of such tissues. However, they do not give the representation of the heating of the magnetic fluid by itself.

Conclusions
A novel simulation model to simulate the heating of a magnetic fluid in an alternating magnetic field in a laboratory setting, considering the coupling between steady state magnetic and transient thermal calculations was built, and the quality of the model was proved with the presented tests. The main issue of the presented model was the appropriate determination of the parameters used for spatial quantities, determination of material's properties, determination of the radiation and convection face regions, and adjustment of equations for coupling between the thermal calculations. All issues were solved. The model was validated with the measurement of a commercial magnetic fluid. The results show a 4.2-24.9% divergence between measurement and simulation in the SAR values.
An analysis of how changing K eff affects the simulation results shows that the accuracy of the calculations depends heavily on the accurate selection of the effective anisotropy constant. Therefore, it is essential how this constant is set. The calculation results show the increase in heating power with rising K eff up until a certain value, where it starts to fall again (Figures 16 and 17). Additional analysis ( Figure 18) shows that the heating power dependency on the K eff reaches a maximum at a certain value (around 30 kJ/m 3 for the magnetic field parameters used in the analysis) at lower temperatures (0-180 • C). With higher K eff values the heating power starts falling, but it starts rising with temperature and reaches a maximum value at around 215 • C.
Magnetic field strength dependence ( Figure 19) shows an increase in time-temperature curves with the increase of the magnetic field strength amplitude. The increase in heating power is more evident in Figure 20, where the initial SAR dependence on magnetic field strength amplitude is presented.
Frequency dependence ( Figure 21) shows an increase in time-temperature with increasing frequency. The heating power increase with the frequency is shown in Figure 22.
Such a model can be used for assessing the heating performance of a magnetic fluid with selected parameters. It can also be used to search for the optimal parameters required to design an optimal magnetic fluid. A model like this has several advantages. It is simple, many simulations can be carried out in a relatively short amount of time, comparisons between fluids with different parameters can be made, and it gives the intrinsic heating properties of the magnetic fluid.
In future work, our research will focus on the development and use of another mathematical model for simulating the magnetic fluid heating process, or simulating another material heating using a similar approach and carrying out the corresponding measurements.