# Comparison of 3D Adaptive Remeshing Strategies for Finite Element Simulations of Electromagnetic Heating of Gold Nanoparticles

Authors: Fadhil Mezghani; Dominique Barchiesi; Abel Cherouat; Thomas Grosges; Houman Borouchaki
Journal: Advances in Mathematical Physics
 (2015)

---

## Abstract

The optical properties of metallic nanoparticles are well known, but the study of their thermal behavior is in its infancy. However the local heating of surrounding medium, induced by illuminated nanostructures, opens the way to new sensors and devices. Consequently the accurate calculation of the electromagnetically induced heating of nanostructures is of interest. The proposed multiphysics problem cannot be directly solved with the classical refinement method of Comsol Multiphysics and a 3D adaptive remeshing process based on ana posteriori error estimator is used. In this paper the efficiency of three remeshing strategies for solving the multiphysics problem is compared. The first strategy uses independent remeshing for each physical quantity to reach a given accuracy. The second strategy only controls the accuracy on temperature. The third strategy uses a linear combination of the two normalized targets (the electric field intensity and the temperature). The analysis of the performance of each strategy is based on the convergence of the remeshing process in terms of number of elements. The efficiency of each strategy is also characterized by the number of computation iterations, the number of elements, the CPU time, and the RAM required to achieve a given target accuracy.


---

## Body

## 1. Introduction

Noble metals nanoparticles, particularly gold nanoparticles, are known as efficient light absorbers [1–3]. Their fascinating optical properties have been studied intensively. Recently, the heat production by these nanoparticles under optical illumination has also aroused much interest [4]. Indeed, an optically excited nanoparticle absorbs a part of the electromagnetic energy that can heat not only the particle itself but also its surrounding [5]. This heating effect becomes particularly strong if the incident electromagnetic wave frequency is near the plasmon resonance frequency of the nanoparticle [6, 7]. The increase of temperature in gold nanoparticles has a variety of applications in nanotechnology, biology, chemistry, and medicine (such as the photothermal cancer therapy using thermal necrosis of tumor cells) [4, 8] and the drug delivery with the remote release of drugs from a capsule containing gold nanoparticles when excited by a laser source [9–11].Owing to these applications, accurate numerical studies of the heating of gold nanoparticle induced by a laser source become pertinent. The problem consists in solving two partial differential equations (PDEs) corresponding to both physics, the result of the electromagnetic problem being the source of heat.The finite element method is widely used for the solution of the partial differential equations, particularly the Helmholtz equation [12] and the heat equation [13]. The mesh quality in the finite element method allows the control of the accuracy of the numerical solutions and the computational memory [14, 15]. The adaption is a powerful technique that can produce an optimal mesh for the finite element calculation [15, 16]. The remeshing methods that produce a new mesh at each iteration of the remeshing loop distinguish themselves from the refinement methods which overall maintain the old mesh but either move some amount of nodes (r-method), increase the degree of interpolation polynomials (p-method), or increase the mesh density by inserting nodes (h-method) or perform a combination of the mentioned methods. The objective of the mesh coarsening is to meet the accuracy requirements fixed by an error estimator that can be physical or geometrical [17, 18]. In this paper, a h-adaptive remeshing process uses a physical error estimator based on ana posteriori p1-interpolation of the physical solution to increase the accuracy of the solution while ensuring the convergence of calculation. Such ana posteriori estimator uses the interpolation of the physical solution computed at the previous step to construct a new mesh that respects the geometry of objects defined in the physical problem. The optimization of the mesh is governed by the Hessian of the interpolation of the solution at the previous step of remeshing. Let us note that the classical h-method inserting nodes on elements edges and therefore dividing each cell of the initial mesh, keeping constant the position of nodes at the previous step of refinement, is not able to solve the investigated problem.Three strategies governing the remeshing are proposed. Each strategy differs from the others by the physical solution that is used for the remeshing process. The accuracy is targeted either on both the electric field and the temperature, on the temperature, or on a linear combination of both, leading to different behavior of the algorithms.Section2 presents the description of the photothermal model and its formulation for finite element method. Section 3 is devoted to both the remeshing adaptive method and the strategies of remeshing using three different choices of physical solutions. In Section 4 numerical results obtained from the three strategies of remeshing are discussed and compared before concluding.


## 2. Physics: Electromagnetic and Thermic Problems

In this section we consider a gold nanoparticle of relative permittivityϵ1 (complex number) and thermal conductivity λ1. The gold nanoparticle is immersed in a uniform dielectric medium with a relative permittivity ϵ2 and a thermal conductivity λ2 (Figure 1). Both media are considered linear, homogeneous, nonmagnetic, and isotropic. The external medium is supposed to be nonabsorbing; therefore the imaginary part of the relative permittivity vanishes: I(ϵ2)=0.Figure 1
Metallic nanosphere immersed in air and illuminated by a monochromatic linearly polarized plane wave.The thermal conductivity of the metallic gold nanoparticle depends on the diameter2R, on the mean free path of electrons in the solid λe,b, and on the thermal conductivity of the bulk metal λbulk. It is can be evaluated from the following relationship [19]:(1)λ1=λbulk1+λe,b/2R.When the radiusR tends toward infinity the thermal conductivity tends toward that of bulk. On the contrary, if 2R is equal to the mean free path of electrons, the thermal conductivity is half of the bulk one. This correction of the thermal properties takes into account the nanometric size of the particles and allows using the heat equation in a wider domain of validity [19].In this study, the spherical nanoparticle of radiusR is illuminated by a monochromatic polarized plane wave Ei=E0exp(ȷωt) with E0 being the amplitude of the incident electric field, ȷ=(-1)1/2 the pure imaginary number of modulus 1 (imaginary unit), and ω=k0c the angular frequency, where k0=2π/λ0 is the modulus of the wave vector of the illumination, λ0 is the wavelength in vacuum, and c is the speed of light in vacuum.In the following, we present the partial derivative equations associated with the multiphysics coupled problem. The interaction between light and matter is governed by the Helmholtz equation derived from the Maxwell equations [12]. The thermic problem is governed by the stationary heat equation [13]. The term of coupling between the electromagnetic and thermal problems is proportional to the square of the electric field inside the metal E2. Let us consider both media successively:(i)
Inside the particle, the electric field E1 satisfies(2)∇×∇×E1-k02ϵ1E1=0,where(∇×) is the curl operator. The electromagnetic power absorbed by the nanosphere is converted to thermal energy. It is considered as the heat source for the thermic problem. It originates from the Joule effect and the temperature inside the particle T1 is given by [13, 20](3)∇·λ1∇T1+ω2ϵ0Iϵ1E12=0,withϵ0 being the permittivity of vacuum; (∇·) and ∇ are the divergence and gradient operators, respectively.(ii)
Outside the particle, the electric field is solution of(4)∇×∇×E2-k02ϵ2E2=0,and the temperature outsideT2 is solution of(5)∇·λ2∇T2=0.To have a full description of the physical problem in a strong formulation of the finite element problem, two boundary conditions must be imposed:(i)
At the interfaceΓ12 between two materials, the boundary conditions correspond to the continuity of the tangential component of the electric field:(6)n12×E2-E1=0,withn12 being the outgoing normal vector. The temperature is supposed to be continuous crossing this boundary T1=T2 and also the heat flux:(7)∂∂n12λ1∇T1=∂∂n12λ2∇T2.(ii)
At the external borderΓ0 of the calculation domain Ω, provided the radiation corresponds to the free propagation of the scattered field and plane wave illumination [12],(8)∇×E2=ȷk0n12×E2-Ei+k0×Eiand a Dirichlet condition for temperature isT0=293.15K (ambient temperature).At the external borderΓ0 of the calculation domain Ω, the imposed conditions are approximations that are relevant if the size of Ω is large enough. For each investigated problem a preliminary numerical study enables finding the minimum size of the domain that leads to stability of results.


## 3. The Adaptive Remeshing and Strategies

In the following, the 3D adaptive remeshing method and the three strategies of remeshing are described.


### 3.1. The Adaptive Remeshing Process

The finite element method is used to solve multiphysics problems in complex geometries [13]. It is based on the subdivision of the computational domain into elements, constituting the mesh. The accuracy of the solution depends on both the shape and the quality of the meshing of the domain. To reach a given accuracy of the solution, a remeshing loop increases the number of elements in regions of Ω where the physical solution exhibits high gradients. At the same time, the density of elements decreases in regions of slight variations of the solution. Nevertheless the total number of elements is always increasing from a step of remeshing to the next one, the minimum size of elements being fixed to capture the variations of the solution. For example, in electromagnetism, the minimum size of elements must be smaller than a fraction of the wavelength. However the calculation on a heavier meshing requires more memory and more computational time. Hence, the number of remeshing steps required to reach a given accuracy is a critical parameter. Therefore the stop criterion for the remeshing loop is the stability of the number of elements: the calculation stops when the relative evolution δ of the number of elements Ni from a step i to the next one is not greater than 10%:(9)δ=Ni+1-NiNi<0.1.This criterion requires the convergence of the remeshing process.The first step of the remeshing method consists in generating a coarse meshing for the geometry of the domain. After a first calculation of the solution, an error estimator is used to generate the physical mapMi [15, 17]. The role of the error estimator is to indicate if the mesh density is too fine or too coarse in some areas of Ω. In this paper we use ana posteriori error estimator based on the Hessian of the interpolation of the physical solution computed at the previous step i-1 of remeshing. The error estimator is deduced from the local error norm of the physical solution. A new mesh of the domain Ω in R3, which respects Mi, is therefore produced. The calculation of the solution and the remeshing process are repeated within a loop until reaching the convergence of the solution towards a fixed precision. The relation between the exact physical solution u, the computed finite element solution uh, and the interpolation field Πhu is as follows [15, 17]: (10)u-uh≤Cu-Πhu,where C is a constant depending on the computational domain. The control of the interpolation error allows the control of the finite element error as shown in the following. Assume a mesh element k with size hk [15, 17]:(11)∀x∈k,ux-Πhux∞≤932Mhk,with M=maxx∈k(max∥ν∥=1|νTHu(x)ν|), where Hu(x) is the Hessian of u at x and ν is any unit direction. Therefore, for controlling the p1-interpolation error as a function of the mesh size, a sufficient condition is(12)hk≤432Mɛ,where ɛ is the maximum error tolerated on the physical solution. Actually this condition implies(13)∀x∈k,ux-Πhux∞≤ɛ.In practice, the exact solution is unknown; therefore the Hessian of the exact solution is approximated by the Hessian of the computed solution uh(x). Therefore the condition on the maximum error becomes(14)∀x∈k,uhx-Πhuhx∞≤ɛ.A solution is acceptable if the local error does not exceed ɛ and if the elements size h respects hmin⁡≤hk≤hmax⁡ in each element of the mesh.The method of remeshing allows generating a new mesh whose density of elements is adapted in areas captured by the error estimator. This formalism allows a better handling of both geometry and topology of the physical solution than the classical refinement. The method is implemented through a modified Delaunay kernel [15, 17].In the following we compare three strategies of meshing based on the different choices of solutions of the multiphysics problem.


### 3.2. Optimization StrategyA1: Controlling Sequentially the Electric Field Intensity and the Temperature

In the model, the electromagnetic and thermic problems are only coupled through the source terms in (3) and therefore both problems can be solved sequentially (weak coupling). The most intuitive method for remeshing is based on the control of accuracy of both problems. Consequently, a first remeshing loop works on the field intensity E2. When target of accuracy is reached for the intensity, the temperature T is calculated in a second remeshing loop. The final mesh of the first loop is used as an initial mesh for the second loop. Algorithm 1 describes optimization strategy A1.Algorithm 1:Algorithm of sequential adaptive remeshing strategy A1.
(1) Given a coarse initial mesh of domainΩ at step i=0: Mi=0(Ω)
(2) loop1:
(3)  Compute‖E‖2 associated with Mi(Ω)
(4)  Define a new discrete metric mapMi on Mi(Ω)
with a posteriori error estimator from ‖E‖2
(5)  Remeshing of the domain with respect to
the metric field (Mi(Ω),Mi)
(6)  ifδ>0.1 (Equation (9))
(7)return to (2) and i=i+1
(8)  else
(9) end loop1i=i+1
(10) loop2:
(11)   ComputeT associated with Mi(Ω)
(12)   Define a new discrete metric mapMi on Mi(Ω)
with a posteriori error estimator from T
(13)   Remeshing of the domain with respect to
the metric field (Mi(Ω),Mi)
(14)   ifδ>0.1 (Equation (9))
(15)  return to (10) and i=i+1
(16)   else
(17) end loop2The electric field near metallic surfaceΓ12 is known to exhibit high gradients [16, 21] unlike the temperature. Consequently the use of the mesh for the electric field is awaited to keep too many elements for the computation of temperature.


### 3.3. StrategyA2: Controlling the Temperature Only

The target being the temperature computation, another intuitive strategy of remeshing considers only this unknown as target of accuracy for remeshing. The optimization loopA2 is described in Algorithm 2.Algorithm 2:Algorithm of adaptive remeshing strategy A2.
(1) Given a coarse initial meshMi=0(Ω)
(2) loop1:
(3)  Compute‖E‖2 and T associated with Mi(Ω)
(4)  Define a new discrete metric mapMi on Mi(Ω)
with a posteriori error estimator from T
(5)  Remeshing of the domain with respect to the metric field(Mi(Ω),Mi)
(6)  ifδ>0.1 (Equation (9))
(7)return to (2) and i=i+1
(8)  else
(9) end loop1The result may be less accurate but is expected to require a smaller number of elements and iterations.


### 3.4. Optimization StrategyA3: Controlling a Linear Combination of Normalized Intensity and Temperature

To generate adequate mesh for both electric fieldE2 and temperature T, the normalized electric field En2 and the normalized temperature Tn are used to construct a linear combination Fn of normalized field as follows:(15)Fn=En2+Tn2,where(16)En2=E2-maxΩ⁡E2maxΩ⁡E2-minΩ⁡E2,Tn=T-maxΩ⁡TmaxΩ⁡T-minΩ⁡T.The size map relative of this new field contains both the areas of strong variation of the electric field and temperature. The optimization loopA3 is illustrated in Algorithm 3.Algorithm 3:Algorithm of adaptive remeshing strategy A3.
(1) Given a coarse initial meshMi=0(Ω)
(2) loop1:
(3)  Compute‖E‖2 and T associated with Mi(Ω)
(4)  ComputeFn (Equation (15))
(5)  Define a new discrete metric mapMi on Mi(Ω)
with a posteriori error estimator from Fn
(6)  Remeshing of the domain with respect to
the metric field (Mi(Ω),Mi)
(7)  ifδ>0.1 (Equation (9))
(8)return to (2) and i=i+1
(9)  else
(10) end loop1The expected result of strategyA3 is a decrease of remeshing steps required to reach accuracy, by taking into account topological properties of both electric field and temperature. The balance between errors between electromagnetic and thermic solution through Fn could lead to a convergence better than that of A2 and to spare computational time compared to A1.


## 3.1. The Adaptive Remeshing Process

The finite element method is used to solve multiphysics problems in complex geometries [13]. It is based on the subdivision of the computational domain into elements, constituting the mesh. The accuracy of the solution depends on both the shape and the quality of the meshing of the domain. To reach a given accuracy of the solution, a remeshing loop increases the number of elements in regions of Ω where the physical solution exhibits high gradients. At the same time, the density of elements decreases in regions of slight variations of the solution. Nevertheless the total number of elements is always increasing from a step of remeshing to the next one, the minimum size of elements being fixed to capture the variations of the solution. For example, in electromagnetism, the minimum size of elements must be smaller than a fraction of the wavelength. However the calculation on a heavier meshing requires more memory and more computational time. Hence, the number of remeshing steps required to reach a given accuracy is a critical parameter. Therefore the stop criterion for the remeshing loop is the stability of the number of elements: the calculation stops when the relative evolution δ of the number of elements Ni from a step i to the next one is not greater than 10%:(9)δ=Ni+1-NiNi<0.1.This criterion requires the convergence of the remeshing process.The first step of the remeshing method consists in generating a coarse meshing for the geometry of the domain. After a first calculation of the solution, an error estimator is used to generate the physical mapMi [15, 17]. The role of the error estimator is to indicate if the mesh density is too fine or too coarse in some areas of Ω. In this paper we use ana posteriori error estimator based on the Hessian of the interpolation of the physical solution computed at the previous step i-1 of remeshing. The error estimator is deduced from the local error norm of the physical solution. A new mesh of the domain Ω in R3, which respects Mi, is therefore produced. The calculation of the solution and the remeshing process are repeated within a loop until reaching the convergence of the solution towards a fixed precision. The relation between the exact physical solution u, the computed finite element solution uh, and the interpolation field Πhu is as follows [15, 17]: (10)u-uh≤Cu-Πhu,where C is a constant depending on the computational domain. The control of the interpolation error allows the control of the finite element error as shown in the following. Assume a mesh element k with size hk [15, 17]:(11)∀x∈k,ux-Πhux∞≤932Mhk,with M=maxx∈k(max∥ν∥=1|νTHu(x)ν|), where Hu(x) is the Hessian of u at x and ν is any unit direction. Therefore, for controlling the p1-interpolation error as a function of the mesh size, a sufficient condition is(12)hk≤432Mɛ,where ɛ is the maximum error tolerated on the physical solution. Actually this condition implies(13)∀x∈k,ux-Πhux∞≤ɛ.In practice, the exact solution is unknown; therefore the Hessian of the exact solution is approximated by the Hessian of the computed solution uh(x). Therefore the condition on the maximum error becomes(14)∀x∈k,uhx-Πhuhx∞≤ɛ.A solution is acceptable if the local error does not exceed ɛ and if the elements size h respects hmin⁡≤hk≤hmax⁡ in each element of the mesh.The method of remeshing allows generating a new mesh whose density of elements is adapted in areas captured by the error estimator. This formalism allows a better handling of both geometry and topology of the physical solution than the classical refinement. The method is implemented through a modified Delaunay kernel [15, 17].In the following we compare three strategies of meshing based on the different choices of solutions of the multiphysics problem.


## 3.2. Optimization StrategyA1: Controlling Sequentially the Electric Field Intensity and the Temperature

In the model, the electromagnetic and thermic problems are only coupled through the source terms in (3) and therefore both problems can be solved sequentially (weak coupling). The most intuitive method for remeshing is based on the control of accuracy of both problems. Consequently, a first remeshing loop works on the field intensity E2. When target of accuracy is reached for the intensity, the temperature T is calculated in a second remeshing loop. The final mesh of the first loop is used as an initial mesh for the second loop. Algorithm 1 describes optimization strategy A1.Algorithm 1:Algorithm of sequential adaptive remeshing strategy A1.
(1) Given a coarse initial mesh of domainΩ at step i=0: Mi=0(Ω)
(2) loop1:
(3)  Compute‖E‖2 associated with Mi(Ω)
(4)  Define a new discrete metric mapMi on Mi(Ω)
with a posteriori error estimator from ‖E‖2
(5)  Remeshing of the domain with respect to
the metric field (Mi(Ω),Mi)
(6)  ifδ>0.1 (Equation (9))
(7)return to (2) and i=i+1
(8)  else
(9) end loop1i=i+1
(10) loop2:
(11)   ComputeT associated with Mi(Ω)
(12)   Define a new discrete metric mapMi on Mi(Ω)
with a posteriori error estimator from T
(13)   Remeshing of the domain with respect to
the metric field (Mi(Ω),Mi)
(14)   ifδ>0.1 (Equation (9))
(15)  return to (10) and i=i+1
(16)   else
(17) end loop2The electric field near metallic surfaceΓ12 is known to exhibit high gradients [16, 21] unlike the temperature. Consequently the use of the mesh for the electric field is awaited to keep too many elements for the computation of temperature.


## 3.3. StrategyA2: Controlling the Temperature Only

The target being the temperature computation, another intuitive strategy of remeshing considers only this unknown as target of accuracy for remeshing. The optimization loopA2 is described in Algorithm 2.Algorithm 2:Algorithm of adaptive remeshing strategy A2.
(1) Given a coarse initial meshMi=0(Ω)
(2) loop1:
(3)  Compute‖E‖2 and T associated with Mi(Ω)
(4)  Define a new discrete metric mapMi on Mi(Ω)
with a posteriori error estimator from T
(5)  Remeshing of the domain with respect to the metric field(Mi(Ω),Mi)
(6)  ifδ>0.1 (Equation (9))
(7)return to (2) and i=i+1
(8)  else
(9) end loop1The result may be less accurate but is expected to require a smaller number of elements and iterations.


## 3.4. Optimization StrategyA3: Controlling a Linear Combination of Normalized Intensity and Temperature

To generate adequate mesh for both electric fieldE2 and temperature T, the normalized electric field En2 and the normalized temperature Tn are used to construct a linear combination Fn of normalized field as follows:(15)Fn=En2+Tn2,where(16)En2=E2-maxΩ⁡E2maxΩ⁡E2-minΩ⁡E2,Tn=T-maxΩ⁡TmaxΩ⁡T-minΩ⁡T.The size map relative of this new field contains both the areas of strong variation of the electric field and temperature. The optimization loopA3 is illustrated in Algorithm 3.Algorithm 3:Algorithm of adaptive remeshing strategy A3.
(1) Given a coarse initial meshMi=0(Ω)
(2) loop1:
(3)  Compute‖E‖2 and T associated with Mi(Ω)
(4)  ComputeFn (Equation (15))
(5)  Define a new discrete metric mapMi on Mi(Ω)
with a posteriori error estimator from Fn
(6)  Remeshing of the domain with respect to
the metric field (Mi(Ω),Mi)
(7)  ifδ>0.1 (Equation (9))
(8)return to (2) and i=i+1
(9)  else
(10) end loop1The expected result of strategyA3 is a decrease of remeshing steps required to reach accuracy, by taking into account topological properties of both electric field and temperature. The balance between errors between electromagnetic and thermic solution through Fn could lead to a convergence better than that of A2 and to spare computational time compared to A1.


## 4. Numerical Results and Discussion

The 3D adaptive remeshing process is applied to compute the heating in the nanosphere and in its immediate vicinity. The gold nanoparticle has complex number relative permittivityϵ1=-11+1.33ȷ (λ0=632.8 nm) as in a previously published 2D study [13]. The value of the relative permittivity of metals is subject to high variability depending on the source of data. It can be computed for each wavelength in the visible domain from [22]: for Johnson and Christy's data it would be ϵ1=-11.29+1.67ȷ and from Palik's handbook ϵ1=-10.61+1.34ȷ. The imaginary part of the relative permittivity obtained from both references is about the same and governs mainly the source of heating. The real parts also are close to each other and therefore the following results must be interpreted in terms of order of magnitude. The thermal conductivity of air is λ2=0.026Wm-1K-1 and that of gold nanoparticle is λ1=131.33Wm-1K-1. The illuminating field is E0=2.667.106(1,0,0)V/m along the x-axis.The radius of the spherical nanoparticle isR=20 nm. The radius of the domain Ω of calculation is RΩ=1400 nm ensuring a stability of electric field and temperature calculation. This problem cannot be solved with the classical refinement method of Comsol Multiphysics. It leads to memory overflow since step 3 of refinement. Both the high ratio RΩ/R and the difference between the topological properties of the electromagnetic field and of the temperature produce an overflow of RAM, from the third iteration of refinement, unlike the remeshing method. We investigate the convergence of the remeshing process, the computational time, and the mesh properties as a function of the above described strategies.


### 4.1. Computational Convergence

To compare the convergence of the remeshing, the parameterhmin is fixed to 0.35 for all strategies. The maximum tolerance on error of the physical solution ɛ is gradually decreased as a function of the remeshing step until reaching 0.005 at step 4: ɛ∈{0.1;0.05,0.01,0.005}. The objective of the calculations is the accurate determination of the temperature; Figures 2(a)-2(b) show the convergence of the strategies. The temperature converges rapidly (after three steps of remeshing) with a maximum shift of 1% between all strategies (Figure 2(a)). Strategy A2 exhibits the best stability, and A3 reaches the convergence criterion at remeshing step i=5. The temperature computed from the initial geometric mesh (i=0) is only 2% overestimated, but a computational result is reliable only if a convergence study is made. The different behavior of strategies is mainly due to the electric field topology exhibiting high gradients. Figure 2(b) illustrates the best convergence of A3 while A2 reaches the stability only after 6 steps of remeshing. Strategy A1 leads to a convergence tray since step 5 but seems to diverge at step 8. The large computing time requires a manual stop of the remeshing to this step.Figure 2
Evolution of the relative maximum of the temperaturemaxΩ(T-T0) (a) and of the intensity of electric field maxΩE2/E02 (b), using the adaptive remeshing method and strategies A1, A2, and A3.
(a)
(b)The maximum overΩ of each solution is a very strong constraint for the convergence of numerical schemes. However the quality of the restoration of topologies of solutions is also of interest. As an illustration, Figure 3 shows a cut of electric field relative intensity E2/E02 and of the temperature elevation ΔT=T-T0 obtained from strategy A2 that seems to be the fastest one. The topology of the relative intensity is characterized by high gradient near the surface of the particle around the direction x of the incoming field (Figure 3) [16]. On the contrary the topology of temperature is smoother (Figure 3). The particle is clearly an isotropic source of heating for the surrounding medium and the continuity of temperature prevents high gradients. The figures obtained from strategies A1 and A3 are indistinguishable from the latter: the relative error between the maps being smaller than 1.5% over Ω.Figure 3
Cut in the median planey=0 of the distribution of the intensity of the electrical field E2/E02: (a) of the temperature T and (b) after adaption with A2.
(a)
(b)At this stage strategyA3 seems to be the best numerical compromise for the accuracy of both the electric field and the temperature. However more specific studies of the computational time and RAM requirements are required before concluding.


### 4.2. Computational Time

The numerical stability of strategies can be evaluated from the number of required steps to reach the convergence criterionδ<10% (9). The number of iteration steps required to reach this criterion is the lowest for A3 and therefore this method could be preferred to A1 that does not converge at the end. Indeed, considering this criterion, A3 is the most efficient method. Figure 4 also reveals the divergence of A1 at iteration 8. Using Fn (15) in A3 balances the influence of both topologies of solutions.Figure 4
Evolution of the computational time for each step with strategiesA1, A2, and A3. The cumulate computational time is about 11 hours for A1, 2 hours for A2, and 3 hours for A3.On the other side, if the computation of temperature must be repeated for the optimization of structures [23] or in an inverse problem [24], the computational time is critical. The cumulative time of computation is linked to the number of steps to reach the convergence criterion. Indeed it is typically about 11 hours for A1, 2 hours for A2, and 3 hours for A3. Consequently, in case of massive use of the multiphysics model, strategy A2 could be the best candidate, even if the accuracy on the temperature is slightly relaxed.The computational time includes the duration of remeshing but is dominated by the inversion time that leads to the solutions. However the RAM used for calculation is a global and limitative criterion of feasibility of the computation.


### 4.3. Memory Usage

An important characteristic of a numerical method is the global memory usage which depends on the number of degrees of freedom and on the number of nodes. The number of degrees of freedom induces the size of the linear system that has to be solved in finite element method. At each step, the optimal meshing is constructed automatically with the adaptive remeshing process, by using the same solver to evaluate the solution. Figures5(a)-5(b) show the difference between the three strategies. Strategy A2 requires the smallest RAM requirement for both electromagnetic and thermic problems. Typically, at the end of the execution of the program, strategy A2 requires 40% less RAM. This property confirms the previous result about the computational time: a smallest linear system is solved faster than a biggest one, even if the number of steps before convergence is one more.Figure 5
Evolution of the RAM for the computation of the electric fieldE (a) and of the temperature T (b) with strategies A1, A2, and A3.
(a)
(b)According to the RAM usage, strategyA2 is the most efficient, but the mesh properties that are discussed in the next subsection could give another point of view.


### 4.4. Mesh Properties

Unnecessary remeshing can be avoided for a given accuracy of computation. This property is discussed in the following. First the convergence of the method, based on the limited increase of number of elementsδ, is shown in Figure 6.Figure 6
Evolution of the number of elements with strategiesA1, A2, and A3.The evolution of the number of elements differs according to the strategy. StrategyA1 does not stop before step 8. This behavior is due to the unnecessary increase of number of elements to describe the gradient of the electric field, whereas the temperature is about the same since step i=3 (Figure 2). On the contrary strategies A2 and A3 stop at steps i=6 and i=5, respectively. Strategy A2 requires twice fewer elements than the others. Strategy A3 is a little bit more efficient than A1. However the accurate plot of solutions in the whole domain Ω of computation is another aspect of the efficiency of the remeshing method. Therefore cuts of the solution and isosurface (that are very sensitive to small variations of solution) give information on this property.Figures7(a)-7(b) show a cut of the initial mesh based only on geometrical properties of objects. This first mesh is the starting point of each strategy. The high density of meshing in the vicinity of the x-direction in Figure 7(c) shows that the remeshing process is able to capture gradients in the electromagnetic problem [16]. The mesh appears to be slightly modified in Figure 7(d) except in the vertical direction showing a tendency to isotropy, a property of temperature. On the contrary this tendency can be clearly seen in Figure 7(e). Moreover, the target being the temperature accuracy, the meshing is much less dense near the surface of the sphere and the topological properties of the intensity do not appear explicitly in the mesh. The density of the elements using A3 strategy reveals equally the handling of the topological properties of both the electric field and the temperature (Figure 7(f)).Figure 7
Cut in the median planey=0: initial mesh (a), 3D view of the mesh (b), and final mesh adapted on E2 with A1(Mi=8(Ω)) (c), adapted on T with A1(Mi=8(Ω)) (d), adapted on T with A2(Mi=6(Ω)) (e), and adapted on Fn with A3(Mi=5(Ω)) (f).
(a)
(b)
(c)
(d)
(e)
(f)The adequate meshing density in the whole domain of calculation can be geometrically characterized by the isosurface of the solution. The isosurfaces are known to be highly sensitive to the discretization of space. Figure8 depicts the isosurfaces of the electric field intensity E2 and of the temperature T associated with the initial mesh Mi=0(Ω) (Figure 7(a)). The bad quality of the representation of solutions E2 (Figure 8(a)) and T (Figure 8(b)) clearly appears. Therefore the initial mesh, based on the geometry of objects, is not sufficient.Figure 8
Isosurfaces of the intensity of electric field (a) and temperature (b) associated with the initial mesh (N0=289605 elements).
(a)
(b)The comparison between isosurfaces obtained from the three strategies at the end of the loops illustrates both the relative accuracy of each method and their ability to get smooth solutions and therefore a good description of both the intensity and the temperature overΩ. Figures 9, 10, and 11 show that the isosurfaces for temperature are almost indistinguishable unlike the isosurfaces for intensity. The electric field intensities obtained from strategies A1 and A3 are close to each other even if the number of elements is smaller for A3. Even if the electric field isosurfaces obtained from strategy A2 appear to be less accurate, the isosurfaces are localized in the same region of Ω. These results show that if the target is only the temperature elevation, A2 is the cheapest strategy.Figure 9
Isosurfaces of the intensity of electric field (a) and temperature (b), after adaptive remeshing withA1 strategy (N8=3380792 tetrahedra).
(a)
(b)Figure 10
Isosurfaces of the intensity of electric field (a) and temperature (b), after adaptive remeshing withA2 strategy (N6=1097595 elements).
(a)
(b)Figure 11
Isosurfaces of the intensity of electric field (a) and temperature (b), after adaptive remeshing withA3 strategy (N5=2586902 elements).
(a)
(b)Table1 summarizes the different properties of each strategy.Table 1
Summary of the efficiency of each strategy,A1, A2, and A3, as a function of the investigated criterion.
A
1
A
2
A
3
RAM
− −
+ +
−
Computational time
− −
+ +
+
Accuracy on electric field
+
−
+ +
Accuracy on temperature
+ +
−
+ +
Convergence (δ criterion)
− −
+
+ +


## 4.1. Computational Convergence

To compare the convergence of the remeshing, the parameterhmin is fixed to 0.35 for all strategies. The maximum tolerance on error of the physical solution ɛ is gradually decreased as a function of the remeshing step until reaching 0.005 at step 4: ɛ∈{0.1;0.05,0.01,0.005}. The objective of the calculations is the accurate determination of the temperature; Figures 2(a)-2(b) show the convergence of the strategies. The temperature converges rapidly (after three steps of remeshing) with a maximum shift of 1% between all strategies (Figure 2(a)). Strategy A2 exhibits the best stability, and A3 reaches the convergence criterion at remeshing step i=5. The temperature computed from the initial geometric mesh (i=0) is only 2% overestimated, but a computational result is reliable only if a convergence study is made. The different behavior of strategies is mainly due to the electric field topology exhibiting high gradients. Figure 2(b) illustrates the best convergence of A3 while A2 reaches the stability only after 6 steps of remeshing. Strategy A1 leads to a convergence tray since step 5 but seems to diverge at step 8. The large computing time requires a manual stop of the remeshing to this step.Figure 2
Evolution of the relative maximum of the temperaturemaxΩ(T-T0) (a) and of the intensity of electric field maxΩE2/E02 (b), using the adaptive remeshing method and strategies A1, A2, and A3.
(a)
(b)The maximum overΩ of each solution is a very strong constraint for the convergence of numerical schemes. However the quality of the restoration of topologies of solutions is also of interest. As an illustration, Figure 3 shows a cut of electric field relative intensity E2/E02 and of the temperature elevation ΔT=T-T0 obtained from strategy A2 that seems to be the fastest one. The topology of the relative intensity is characterized by high gradient near the surface of the particle around the direction x of the incoming field (Figure 3) [16]. On the contrary the topology of temperature is smoother (Figure 3). The particle is clearly an isotropic source of heating for the surrounding medium and the continuity of temperature prevents high gradients. The figures obtained from strategies A1 and A3 are indistinguishable from the latter: the relative error between the maps being smaller than 1.5% over Ω.Figure 3
Cut in the median planey=0 of the distribution of the intensity of the electrical field E2/E02: (a) of the temperature T and (b) after adaption with A2.
(a)
(b)At this stage strategyA3 seems to be the best numerical compromise for the accuracy of both the electric field and the temperature. However more specific studies of the computational time and RAM requirements are required before concluding.


## 4.2. Computational Time

The numerical stability of strategies can be evaluated from the number of required steps to reach the convergence criterionδ<10% (9). The number of iteration steps required to reach this criterion is the lowest for A3 and therefore this method could be preferred to A1 that does not converge at the end. Indeed, considering this criterion, A3 is the most efficient method. Figure 4 also reveals the divergence of A1 at iteration 8. Using Fn (15) in A3 balances the influence of both topologies of solutions.Figure 4
Evolution of the computational time for each step with strategiesA1, A2, and A3. The cumulate computational time is about 11 hours for A1, 2 hours for A2, and 3 hours for A3.On the other side, if the computation of temperature must be repeated for the optimization of structures [23] or in an inverse problem [24], the computational time is critical. The cumulative time of computation is linked to the number of steps to reach the convergence criterion. Indeed it is typically about 11 hours for A1, 2 hours for A2, and 3 hours for A3. Consequently, in case of massive use of the multiphysics model, strategy A2 could be the best candidate, even if the accuracy on the temperature is slightly relaxed.The computational time includes the duration of remeshing but is dominated by the inversion time that leads to the solutions. However the RAM used for calculation is a global and limitative criterion of feasibility of the computation.


## 4.3. Memory Usage

An important characteristic of a numerical method is the global memory usage which depends on the number of degrees of freedom and on the number of nodes. The number of degrees of freedom induces the size of the linear system that has to be solved in finite element method. At each step, the optimal meshing is constructed automatically with the adaptive remeshing process, by using the same solver to evaluate the solution. Figures5(a)-5(b) show the difference between the three strategies. Strategy A2 requires the smallest RAM requirement for both electromagnetic and thermic problems. Typically, at the end of the execution of the program, strategy A2 requires 40% less RAM. This property confirms the previous result about the computational time: a smallest linear system is solved faster than a biggest one, even if the number of steps before convergence is one more.Figure 5
Evolution of the RAM for the computation of the electric fieldE (a) and of the temperature T (b) with strategies A1, A2, and A3.
(a)
(b)According to the RAM usage, strategyA2 is the most efficient, but the mesh properties that are discussed in the next subsection could give another point of view.


## 4.4. Mesh Properties

Unnecessary remeshing can be avoided for a given accuracy of computation. This property is discussed in the following. First the convergence of the method, based on the limited increase of number of elementsδ, is shown in Figure 6.Figure 6
Evolution of the number of elements with strategiesA1, A2, and A3.The evolution of the number of elements differs according to the strategy. StrategyA1 does not stop before step 8. This behavior is due to the unnecessary increase of number of elements to describe the gradient of the electric field, whereas the temperature is about the same since step i=3 (Figure 2). On the contrary strategies A2 and A3 stop at steps i=6 and i=5, respectively. Strategy A2 requires twice fewer elements than the others. Strategy A3 is a little bit more efficient than A1. However the accurate plot of solutions in the whole domain Ω of computation is another aspect of the efficiency of the remeshing method. Therefore cuts of the solution and isosurface (that are very sensitive to small variations of solution) give information on this property.Figures7(a)-7(b) show a cut of the initial mesh based only on geometrical properties of objects. This first mesh is the starting point of each strategy. The high density of meshing in the vicinity of the x-direction in Figure 7(c) shows that the remeshing process is able to capture gradients in the electromagnetic problem [16]. The mesh appears to be slightly modified in Figure 7(d) except in the vertical direction showing a tendency to isotropy, a property of temperature. On the contrary this tendency can be clearly seen in Figure 7(e). Moreover, the target being the temperature accuracy, the meshing is much less dense near the surface of the sphere and the topological properties of the intensity do not appear explicitly in the mesh. The density of the elements using A3 strategy reveals equally the handling of the topological properties of both the electric field and the temperature (Figure 7(f)).Figure 7
Cut in the median planey=0: initial mesh (a), 3D view of the mesh (b), and final mesh adapted on E2 with A1(Mi=8(Ω)) (c), adapted on T with A1(Mi=8(Ω)) (d), adapted on T with A2(Mi=6(Ω)) (e), and adapted on Fn with A3(Mi=5(Ω)) (f).
(a)
(b)
(c)
(d)
(e)
(f)The adequate meshing density in the whole domain of calculation can be geometrically characterized by the isosurface of the solution. The isosurfaces are known to be highly sensitive to the discretization of space. Figure8 depicts the isosurfaces of the electric field intensity E2 and of the temperature T associated with the initial mesh Mi=0(Ω) (Figure 7(a)). The bad quality of the representation of solutions E2 (Figure 8(a)) and T (Figure 8(b)) clearly appears. Therefore the initial mesh, based on the geometry of objects, is not sufficient.Figure 8
Isosurfaces of the intensity of electric field (a) and temperature (b) associated with the initial mesh (N0=289605 elements).
(a)
(b)The comparison between isosurfaces obtained from the three strategies at the end of the loops illustrates both the relative accuracy of each method and their ability to get smooth solutions and therefore a good description of both the intensity and the temperature overΩ. Figures 9, 10, and 11 show that the isosurfaces for temperature are almost indistinguishable unlike the isosurfaces for intensity. The electric field intensities obtained from strategies A1 and A3 are close to each other even if the number of elements is smaller for A3. Even if the electric field isosurfaces obtained from strategy A2 appear to be less accurate, the isosurfaces are localized in the same region of Ω. These results show that if the target is only the temperature elevation, A2 is the cheapest strategy.Figure 9
Isosurfaces of the intensity of electric field (a) and temperature (b), after adaptive remeshing withA1 strategy (N8=3380792 tetrahedra).
(a)
(b)Figure 10
Isosurfaces of the intensity of electric field (a) and temperature (b), after adaptive remeshing withA2 strategy (N6=1097595 elements).
(a)
(b)Figure 11
Isosurfaces of the intensity of electric field (a) and temperature (b), after adaptive remeshing withA3 strategy (N5=2586902 elements).
(a)
(b)Table1 summarizes the different properties of each strategy.Table 1
Summary of the efficiency of each strategy,A1, A2, and A3, as a function of the investigated criterion.
A
1
A
2
A
3
RAM
− −
+ +
−
Computational time
− −
+ +
+
Accuracy on electric field
+
−
+ +
Accuracy on temperature
+ +
−
+ +
Convergence (δ criterion)
− −
+
+ +


## 5. Conclusion

A 3D adaptive remeshing process for the generation, adaption, and optimization of meshes has been applied to the field of electromagnetics-thermics coupled problem. The adaptive remeshing process is based on ana posteriori error estimator of p1-interpolation and therefore strongly differs from refinement techniques that do not permit solving the investigated problem. The remeshing method permits changing automatically the element size for minimizing the interpolation error of numerical solutions. However the choice of the physical field used for the remeshing method leads to very different behaviors. In this study we considered E2 and T sequentially (strategy A1), only T (strategy A2), and a linear combination of normalized E2 and T (strategy A3). The three strategies produce accurate solutions but their stability and efficiency differ. The convergence of A1 is not ensured due to the drastic increase of the number of elements required to describe the electric field gradient (but useless for temperature). Strategy A3 is more stable and converges more rapidly than A1 and requires less computational resources. Strategy A2 is much less costly, to the detriment of the accuracy on the temperature and of the ability to plot the electric field. The choice between methods A2 and A3 should therefore be made according to a compromise between accuracy on one hand and constraints on time and memory on the other hand. The nonconvergence (using the criterion of the increase of the number of elements) of classical and natural method A1 should lead to avoiding it for computation of temperature elevation induced by electromagnetic illumination. Such a remeshing process and strategy A3 can be used to investigate the temperature evolution in the vicinity of various nanoparticles used in cancer therapy [25].



---

*Source: 469310-2015-05-24.xml*