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 an a 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.
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 (-method), increase the degree of interpolation polynomials (-method), or increase the mesh density by inserting nodes (-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 -adaptive remeshing process uses a physical error estimator based on an a posteriori p1-interpolation of the physical solution to increase the accuracy of the solution while ensuring the convergence of calculation. Such an a 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 -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.
Section 2 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 (complex number) and thermal conductivity . The gold nanoparticle is immersed in a uniform dielectric medium with a relative permittivity and a thermal conductivity (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: .
The thermal conductivity of the metallic gold nanoparticle depends on the diameter , on the mean free path of electrons in the solid , and on the thermal conductivity of the bulk metal . It is can be evaluated from the following relationship [19]:
When the radius tends toward infinity the thermal conductivity tends toward that of bulk. On the contrary, if 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 radius is illuminated by a monochromatic polarized plane wave with being the amplitude of the incident electric field, the pure imaginary number of modulus (imaginary unit), and the angular frequency, where is the modulus of the wave vector of the illumination, is the wavelength in vacuum, and 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 . Let us consider both media successively:(i) Inside the particle, the electric field satisfies 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 is given by [13, 20] with being the permittivity of vacuum; and are the divergence and gradient operators, respectively.(ii) Outside the particle, the electric field is solution of and the temperature outside is solution of
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 between two materials, the boundary conditions correspond to the continuity of the tangential component of the electric field: with being the outgoing normal vector. The temperature is supposed to be continuous crossing this boundary and also the heat flux:(ii)At the external border of the calculation domain , provided the radiation corresponds to the free propagation of the scattered field and plane wave illumination [12], and a Dirichlet condition for temperature is (ambient temperature).
At the external border 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 from a step to the next one is not greater than :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 map [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 an a posteriori error estimator based on the Hessian of the interpolation of the physical solution computed at the previous step of remeshing. The error estimator is deduced from the local error norm of the physical solution. A new mesh of the domain in , which respects , 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 , the computed finite element solution , and the interpolation field is as follows [15, 17]: where 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 with size [15, 17]:with , where is the Hessian of at and is any unit direction. Therefore, for controlling the p1-interpolation error as a function of the mesh size, a sufficient condition iswhere is the maximum error tolerated on the physical solution. Actually this condition impliesIn practice, the exact solution is unknown; therefore the Hessian of the exact solution is approximated by the Hessian of the computed solution . Therefore the condition on the maximum error becomesA solution is acceptable if the local error does not exceed and if the elements size respects 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 Strategy : 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 . When target of accuracy is reached for the intensity, the temperature 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 .
| ||||||||||||||||||||||||||||||||||||||||||
The electric field near metallic surface 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. Strategy : 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 loop is described in Algorithm 2.
| ||||||||||||||||||||
The result may be less accurate but is expected to require a smaller number of elements and iterations.
3.4. Optimization Strategy : Controlling a Linear Combination of Normalized Intensity and Temperature
To generate adequate mesh for both electric field and temperature , the normalized electric field and the normalized temperature are used to construct a linear combination of normalized field as follows:where
The size map relative of this new field contains both the areas of strong variation of the electric field and temperature. The optimization loop is illustrated in Algorithm 3.
| ||||||||||||||||||||||||
The expected result of strategy 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 could lead to a convergence better than that of and to spare computational time compared to .
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 ( 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 and from Palik's handbook . 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 and that of gold nanoparticle is . The illuminating field is along the -axis.
The radius of the spherical nanoparticle is nm. The radius of the domain of calculation is 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 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 parameter is fixed to for all strategies. The maximum tolerance on error of the physical solution is gradually decreased as a function of the remeshing step until reaching at step 4: . 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 exhibits the best stability, and reaches the convergence criterion at remeshing step . The temperature computed from the initial geometric mesh () 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 while reaches the stability only after 6 steps of remeshing. Strategy 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.
(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 and of the temperature elevation obtained from strategy 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 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 and are indistinguishable from the latter: the relative error between the maps being smaller than 1.5% over .
(a)
(b)
At this stage strategy 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 (9). The number of iteration steps required to reach this criterion is the lowest for and therefore this method could be preferred to that does not converge at the end. Indeed, considering this criterion, is the most efficient method. Figure 4 also reveals the divergence of at iteration 8. Using (15) in balances the influence of both topologies of solutions.
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 , 2 hours for , and 3 hours for . Consequently, in case of massive use of the multiphysics model, strategy 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. Figures 5(a)-5(b) show the difference between the three strategies. Strategy requires the smallest RAM requirement for both electromagnetic and thermic problems. Typically, at the end of the execution of the program, strategy 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.
(a)
(b)
According to the RAM usage, strategy 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.
The evolution of the number of elements differs according to the strategy. Strategy 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 (Figure 2). On the contrary strategies and stop at steps and , respectively. Strategy requires twice fewer elements than the others. Strategy is a little bit more efficient than . 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.
Figures 7(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 -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 strategy reveals equally the handling of the topological properties of both the electric field and the temperature (Figure 7(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. Figure 8 depicts the isosurfaces of the electric field intensity and of the temperature associated with the initial mesh (Figure 7(a)). The bad quality of the representation of solutions (Figure 8(a)) and (Figure 8(b)) clearly appears. Therefore the initial mesh, based on the geometry of objects, is not sufficient.
(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 and are close to each other even if the number of elements is smaller for . Even if the electric field isosurfaces obtained from strategy 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, is the cheapest strategy.
(a)
(b)
(a)
(b)
(a)
(b)
Table 1 summarizes the different properties of each strategy.
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 an a 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 and sequentially (strategy ), only (strategy ), and a linear combination of normalized and (strategy ). The three strategies produce accurate solutions but their stability and efficiency differ. The convergence of 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 is more stable and converges more rapidly than and requires less computational resources. Strategy 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 and 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 should lead to avoiding it for computation of temperature elevation induced by electromagnetic illumination. Such a remeshing process and strategy can be used to investigate the temperature evolution in the vicinity of various nanoparticles used in cancer therapy [25].
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.