Accuracy Assessment of Numerical Dosimetry for the Evaluation of Human Exposure to Electric Vehicle Inductive Charging Systems

In this article, we discuss numerical aspects related to the accuracy and the computational efficiency of numerical dosimetric simulations, performed in the context of human exposure to static inductive charging systems of electric vehicles. Two alternative numerical methods based on electric vector potential and electric scalar potential formulations, respectively, are here considered for the electric field computation in highly detailed anatomical human models. The results obtained by the numerical implementation of both approaches are discussed in terms of compliance assessment with ICNIRP guidelines limits for human exposure to electromagnetic fields. In particular, different strategies for smoothing localized unphysical outliers are compared, including novel techniques based on statistical considerations. The outlier removal is particularly relevant when comparison with basic restrictions is required to define the safety of electromagnetic fields exposure. The analysis demonstrates that it is not possible to derive general conclusions about the most robust method for dosimetric solutions. Nevertheless, the combined use of both formulations, together with the use of an algorithm for outliers removal based on a statistical approach, allows to determine final results to be compared with reference limits with a significant level of reliability.


I. INTRODUCTION
T HE WIDE spreading of wireless power charging technologies has increased the public concern about potential effects of human exposure to electromagnetic fields (EMFs) produced by such equipment. Among them, inductive charging of electric vehicles (EVs) represents a technology in rapid development, being a major step towards a new transport paradigm. Indeed, inductive charging (particularly while driving) has many advantages: one above all, it will make the recharging phase easier and more favorable for the driver. Static inductive power transfer (IPT) is obviously the first stage of this technology, but the market penetration of this technology needs rigorous safety consideration for humans (drivers, passengers, and by-standers) because exposure to magnetic stray fields could be, to some extent, significant in the vicinity of EVs, also considering the levels of electric power involved. This topic is of great relevance nowadays, as demonstrated by the numerous papers recently published in the literature (e.g., [1]- [10]). In order to support the exposure analysis, numerical dosimetry is largely adopted when the estimation of the induced fields in the human body is needed. Indeed, the ICNIRP guidelines [11] require the compliance with basic restrictions (BRs), i.e., induced electric field levels in the human body when the exposure occurs in the presence of highly inhomogeneous magnetic fields. A recent evaluation of the exposure of a human standing around a vehicle charged via IPT, presented in [10], shows that compliance with ICNIRP guidelines in terms of BRs is verified, despite the magnetic flux density exceeds the reference levels.
Tools for numerical dosimetry have had an extraordinary development over the last decade, by virtue of the potentiality of last generation personal computers, the software development (e.g., [12]- [15]), and the availability of highly sophisticated human body models with associated electrical properties of tissues (e.g., [16]- [18]). Nevertheless, the accuracy of the dosimetric results, including the detection of possible numerical artefacts, remains a topic of great interest, particularly when a comparison with exposure limits is needed [13], [19].
The estimation of the reliability of dosimetric simulations is a common task for different exposure scenarios, but the case of inductive charging of electric vehicles is of great relevance because of the significant level of electric power involved at a relatively high frequency (usually around some tens of kilohertz). Moreover, the stray magnetic field is highly inhomogeneous and the exposure potentially involves the entire population, including children and babies.
The study of human exposure related to IPT is quite complex, involving at least two computational challenges. The first one consists in the determination of the stray magnetic field around This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see http://creativecommons.org/licenses/by/4.0/ or inside the EVs, i.e., in an open-boundary domain. In such a volume, the transmitting and receiving coils are close to the region of interest (distance ∼ 1 m, near-field analysis) and thin metallic/magnetic car body structures and related components (e.g., magnetic flux concentrators) significantly affect the local magnetic field distribution. The second computational challenge concerns the evaluation in realistic human models of induced EMFs, whose distribution is a function of the shape of the whole exposed subject. However, the analysis of the entire human body has consequence on the accuracy that can be reached by the simulations, which depends on the resolution of the computational domain and on the number of unknowns, particularly large in highly detailed anatomical models.
The main scope of this article is to investigate the accuracy of dosimetric results in the framework of human exposure to IPT systems, supplementing results already published in the literature for different exposure scenarios [13], [20], [21]. As an example, in [21], six codes are compared based on impedance method, scalar potential finite-difference method, quasi-static finite-difference time domain method, and fastmultipole surface-charge-simulation method for voxel data. After being validated on homogeneous and/or multi-layer spheres, such codes have been compared by analyzing a whole-body model (the Japanese male TARO with 2-mm discretization) exposed to a uniform 0.1 mT magnetic field at 50 Hz. Apart from the impedance method, which provides results significantly higher than the other ones, the results of the other codes exhibit discrepancies up to 25%.
In this article, we focus in particular on the comparison of two further alternative finite element formulations to solve the EMF problem related to the dosimetric study. They are complementary formulations, one based on the electric scalar potential, largely adopted in the literature for dosimetric analysis (see, for example, [12], [22]) and implemented in scientific commercial softwares for dosimetric analysis (see, for example, [23]), and the other based on the electric vector potential, usually applied in the analysis of eddy current problems (e.g., [24]). Attention is focused on the robustness of the electromagnetic solvers based on the two alternative field formulations to the rise of numerical artifacts (unphysical outliers), discussing their merits and drawbacks in terms of accuracy.
In addition to compare the raw results (not filtered) obtained by the two formulations, we investigate the potentialities of an algorithm for the outlier removal, proposed in literature in [20] and here revisited. This algorithm is applied to the raw results obtained by the two considered formulations with the aim of finding, through the intercomparison of the results, a more reliable prediction.
The accuracy assessment of the dosimetric simulations is discussed in terms of the maximum level of induced electric field computed within the body, which is the value that has to comply with the ICNIRP BRs [11]. The related considerations could have some impact on the definition of the most appropriate metric for human exposure assessment in IPT scenarios with respect to the ICNIRP guidelines, as also discussed in [10]. The analysis is developed making reference to a specific IPT system, but most of the considerations here reported go beyond the specific exposure scenario, being applicable to similar low and medium frequency human exposure contexts.
To the authors' knowledge, the intercomparison of complementary finite element formulations, supplemented by the application of algorithms for outliers removal for the determination of the most appropriate metric for human exposure assessment in IPT systems, represents the main novelty of this article with respect to the relevant literature.

A. Model of the Considered IPT System
In the analysis, we make reference to an IPT system suitable for light electric vehicles (transmission power of ∼11 kW at 85 kHz [25]). The system consists of a transmitter coil, drowned under the asphalt, and a receiver coil, placed on-board of the vehicle. Flux concentrators optimize the coupling between primary and secondary circuits. An aluminum sheet structure, usually attached to the receiving coil, mitigates the field inside the vehicle. The presence of the body car is here disregarded. A simplified scheme of the main components is reported in Fig. 1. The transmitter coil has a size of 1.62 m × 0.62 m; the receiver coil has a size of 0.62 m × 0.42 m and it is placed at a distance of 0.20 m. Both coils have ten turns. The aluminum sheet, having size equal to 1.41 m × 0.52 m, is placed at 0.14 m from the receiver coil.

B. Method for the Computation of the Stray Magnetic Field
Taking advantage of the weak-field reaction due to the body's tissues (this approximation introduces an error on the induced electric field lower than 1 ppm, which was estimated using the analytical solution reported in [26]), the analysis of human exposure can be split into two steps.
In the first one, the stray magnetic field H s generated by the IPT system in the body volume is computed, neglecting the presence of the human body but including the effects of the conductive/ferromagnetic components. In the second step, the computed spatial distribution of magnetic field is applied to a highly detailed anatomical model, placed in a given position and discretized into voxels, to evaluate the electric field induced within the body.
The stray magnetic field generated by the transmitting and receiving coils is here computed using the dynamic electromagnetics module of the simulation software Opera-3D [27] able to take into account the effects of flux concentrators and thin metallic structures. Under time-periodic working conditions, the analysis is conveniently performed in the frequency domain. The outputs are the Cartesian components of the magnetic field H s (real and imaginary part) in the spatial positions corresponding to the barycenter of each voxel of the human model.

C. Finite Element Formulations for the Dosimetric Analysis
For the dosimetric simulations, two alternative formulations [12] are here used and compared. In the first formulation, an edge-based finite-element method solver is adopted, assuming  [25]). In the lower figure, a plan view, with the indication of the reference system and the main size (in meters). In both figures, the position of the anatomical body model (here the ViP model Roberta) is also represented, at the distance of 1 m from the center of the reference system. Please note that the receiver coil is centered with respect to the transmitter coil. an electric vector potential T (J = curl T) as the unknown (Tformulation). The weak form equation for the domain Ω (i.e., the human body model) is being σ the spatial-dependent electrical conductivity of tissues, ω the angular frequency of the magnetic field source, μ 0 the air magnetic permeability, j the imaginary unit, and υ the vector test function. Implicit in (1) is the Dirichlet boundary condition n × T = 0, which imposes body-confined induced currents.
Equation (1) is implemented using an ungauged approach, considering all internal edges as unknowns, without introducing tree-cotree separation. This choice increases the number of unknowns but favors the convergence of the iterative solver [28]. Alternatively, a nodal finite-element solver, based on the A−φ formulation, was used assuming the electric scalar po- being w the scalar test function. The Neumann boundary condition grad φ · n = jωA s · n imposes that induced currents are confined within the body. This approach corresponds to the so-called scalar potential finite difference [22]. Also, for this formulation, we have adopted an ungauged implementation, which favors the iterative solver.
In (2), the unperturbed field sources are represented by the magnetic vector potential A s . This later can be directly obtained by Cobham Opera-3D solver or, alternatively, computed by curl-inversion of the H s spatial distribution [29]- [31]. Only the last approach is considered and tested here, because it can be applied also starting from measured maps of the source magnetic field [32].
The dosimetric problems (1) and (2) are solved numerically by using a finite-element discretization into hexahedra coincident with the voxels of the human body model. Due to the large number of unknowns arising from the use of highly detailed anatomical models, the algebraic system of equations is solved through a generalized minimal residual method (GMRES) algorithm [33], which allows avoiding the storage of the FEM matrix. The solution of problems (1) and (2) provides the estimate of the electric field magnitudes induced in the barycenter of each voxel of the anatomical model

III. COMPARISON WITH ANALYTICAL SOLUTION
Before moving to the analysis of anatomical human models, the two formulations have been tested with respect to an analytical solution available for two nested but misaligned cylinders (of infinite length along their axes), having different electrical conductivities [34]. A uniform magnetic field, directed along the cylinder axes (z-axis), produces eddy currents circulating in transversal xy-planes (2D problem). The magnetic effects generated by eddy currents are disregarded, as in the case of biological tissues. The advantage of this configuration with respect to the multilayered spheres, which also provides an analytical solution, is that the spatial distribution of the induced currents determines different interface situations, with prevalence of tangential or normal components of the field quantities, allowing a more severe test for the two complementary formulations considered in this article.
The considered external and internal cylinders have diameters, respectively, equal to 0.4 m and 0.1 m, with 0.1 m misalignment between the centers of the two cylinders. Two situations are simulated: one with the electrical conductivity of the internal cylinder equal to 1 S/m and the external cylinder having conductivity equal to 0.01 S/m (case #1), and the other one with opposite values (case #2). A unitary magnetic flux density is imposed, with a frequency of 85 kHz. In the 3D numerical solutions, where the extent of the cylinders has been limited along the z-direction, the geometry has been meshed with cubic voxels of size equal to 1 mm and 2 mm. The different contrast between the electric conductivities in the two cases determines different spatial distribution of the induced currents, as evidenced in Fig. 2.
Numerical local values of the induced electric field (E num .) are compared with the analytical ones (E anal .) computing the following index errors in the central xy-section of the cylinders: being Δx and Δy the voxel size along the x-and y-axis (transversal plane), respectively. The resulting numerical values are summarized in Table I, together with the spatial localization of the maximum error. For case #1, it results that the maximum deviation in terms of electric field values ranges between 3753 V/m and 4045 V/m (to be compared with the analytical maximum value of 72 000 V/m) for all simulations, with a slightly better performance of the A−φ formulation. This is also confirmed by the L2 error. Also, for case #2, the A−φ formulation performs slightly better, with a maximum deviation in terms of electric field values ranging between 4908 V/m and 7462 V/m (to be compared with the analytical maximum value of 120 000 V/m) for all simulations. It must be noted that the localization of the maximum error is always at the interface between the two cylinders for the A−φ formulation, but it moves to the external boundary for the T formulation in case #1. The localization of the maximum error at the interfaces is due to the different conservation properties implicitly imposed by the two complementary formulations at the interfaces; indeed, while T formulation preserves the normal component of current density J, the A−φ formulation preserves the tangential component of the electric field E. In addition, a geometrical discrepancy occurs in those regions because curved interfaces are unavoidably discretized as staircase surfaces when using voxels. This effect will be discussed in the next Section IV-A.
Finally, it is worth noting that despite the L2 error reduces at the reduction of the voxel size, demonstrating a better global accuracy, the maximum local error slightly increases, because of a reduced averaging effect in smaller voxels. This behavior will be evident also in the successive analysis with the anatomical human models.

A. Numerical/Geometrical Artefacts
Numerical solutions computed on voxel-based models are usually affected by numerical artefacts, as evidenced in the previous paragraph and as discussed in several papers (see, for example, [13], [20], [21]). The main causes of numerical artefacts are the stair-casing error [13], inevitable when using voxelized models, and the contrast between electromagnetic properties of adjacent voxels belonging to different tissues. In addition, local artefacts may occur when the voxelization introduces local contacts not present in the original CAD of the human body; these artefacts appear especially in proximity of the skin (with reference to skin-to-skin contact voxels see the approved Draft of the IEEE C95.1 Standard [35]).
All these effects introduce fictitious outliers, which have to be avoided, particularly when they alter the maximum level of the induced electric field that is the metric to be used for evaluating the compliance with ICNIRP BRs.
To evidence how the stair-casing error is intrinsic in voxelized domains, this effect has been here reproduced in a 2D axisymmetric model problem, which despite its simplicity allows to compare the artefact with respect to a reference solution. In this model problem, a sphere of electrical conductivity equal to 0.1 S/m is embedded within a matrix of conductivity equal to 0.35 S/m. A uniform magnetic field (unitary amplitude) is applied along the z-axis to generate eddy currents in transversal xy-planes. The interface between the two considered tissues is modeled by a smooth (reference solution) and a staircase interface, which reproduces in the 2D model the same discretization of the analogous 3D voxel model. The 2D axisymmetric solutions have been obtained by a conventional finite-element solver [36], based on a triangular mesh, which provides a stable solution in terms of mesh discretization. The 3D solutions to the same problem have been obtained by applying the two formulations described above. The results are compared in Fig. 3, by plotting the magnitude of the E-field along a line in close proximity of the interface between tissues. The local oscillations are well evident in all models with a staircase interface, including the 2D solution with stepped boundary. In other words, the oscillation at the interface is the actual behavior of the electric field in the staircase configuration and the lack of accuracy is essentially due to imperfections of the geometrical model not to inaccuracies of the numerical solutions. It must be noted that, as highlighted in previous studies [14], [20], the use of fine voxel-based meshes (useful to keep high quality in the anatomical description) has the drawback to magnify the presence of outliers.

B. Filtering Techniques
Filtering techniques can be applied on raw E-field (E raw ) datasets to overcome the unavoidable unphysical outliers that affect the 3D solutions based on voxel models. One of these filtering approaches is proposed by ICNIRP [10], which takes the presence of outliers for granted and therefore requires to compare the BRs to the 99th percentile of the induced electric field magnitude (E (99th) ), instead of its absolute maximum. However, previous studies (e.g., [37]- [39]) have put in evidence some faults of such an approach, which may overestimate the presence of outliers in the case of a strongly heterogeneous exposure. For instance, the hand of a person handling an electrical device that produces a magnetic field results to be much more exposed to the field itself than the rest of the body. In this case, the application of the 99th percentile (for the whole body, or for a very extensive specific tissue, like the skin) may introduce an overcorrection and remove genuine "hot-spots," which should be taken into account in the exposure assessment.
To overcome some of the flaws in the metric recommended by ICNIRP [11], another way for removing outliers, more conservative with respect to the 99th percentile, has been proposed in [20] and is here improved to remove some mismatches generated during the extrapolation of the filtered data.
The technique proposed in [20], briefly recalled here for the readers' convenience, requires to sort in ascending order the magnitude of all E-field values, despite their spatial distribution. Typically, the curve of the sorted values is quite smooth, but, in the last part, exhibits a strong upward concavity due to the sudden rise of the highest values (see next Fig. 8). The "gradient" (defined as the difference between two adjacent values divided by their mean) of the last percentile of the sorted values (i.e., the highest values) is then computed. If the frequency distribution of the logarithm of the gradient is almost symmetric and unimodal (i.e., exhibits one local maximum only), the detection point (DP) of the outliers is determined by adding to the mean of the distribution a positive shift, corresponding to three times the standard deviation (for a Gaussian distribution, this choice  covers 99.87 % of the smallest gradient values). An example of the frequency distribution is shown in Fig. 4 considering the E-field values computed by the A−φ formulation for an anatomical model with 1 mm resolution (simulations discussed in the next section). The DP is evidenced in the diagram.
Once the DP has been identified, we search, within the last percentile of the sorted E-field values, the first pair of values that produce a gradient whose logarithm is larger than DP. From these values on, the computed E-field values are considered as outliers and are corrected. In [20], the correction is obtained by fitting the sorted values that are nonoutliers (within the last percentile of all E-field values) using a polynomial function of order 2, and replacing the outliers with the extrapolation of this trend.
However, as shown in Fig. 5, this correction may introduce a significant discontinuity in the sorted values, across the border between the values that do not undergo the correction (identified as nonoutliers) and the values obtained through the polynomial extrapolation. This happens because the interpolating curve does not necessarily pass through the last nonoutliers. The simplest possibility to remove the discontinuity is to subtract the amount of the discontinuity itself from the extrapolated values, so that the latter are forced to match the value of the last nonoutliers. This solution restores the continuity of the curve of the sorted E-field values, but, in general, introduces a sharp change in the slope of such a curve.
In order to avoid all these faults, which determine different values of the maximum filtered electric field, we propose here to calibrate the interpolating curve on the last nonoutliers, whose number is chosen equal to the number of the outliers to be corrected. This choice, which constitutes the main improvement with respect to [20], typically reduces the population of nonoutliers used to identify the polynomial fitting and guarantees a good continuity of the resulting curve, as well as of its first derivative, at the interface between the points kept at their original values and those obtained via extrapolation. In the following, these values are named as E (out) .

A. Models for In-Silico Simulations
We move now to the application of the methods described above to the analysis of human exposure to the IPT system described in Section II-A.
The exposure analysis is performed during the charging phase, when both the transmitter and the receiver coils are energized. In particular, the peak value of the current in the transmitter coil is 50 A, while in the receiver coil, it is 105 A, in quadrature with respect to the primary current. It must be remarked that IPT systems do not work like a standard transformer. Although the efficiency of the up-to-date systems is very high (>90%), the absence of a ferromagnetic core makes the coefficient of coupling low. For this reason, the ratio between the number of turns is not related with the ratio between transmitting and receiving coil current [40].
Different anatomical models, all belonging to the virtual population (ViP) developed by the IT'IS Foundation [16], have been used to check the accuracy of the dosimetric estimation, considering the dielectric properties derived from the IT'IS Database [17]. In particular, we have considered: 1) the 5-years old female "Roberta," with height 1.09 m and weight 17.8 kg, segmented into 65 different tissues; 2) the 34-years male adult "Duke," with height 1.77 m and weight 70.3 kg, segmented into 305 different tissues. The posable version of this anatomical model has been used to obtain, besides the original standing position, also a crouching position. For all models, the electrical conductivity ranges from 0.02 S/m (bones) to 2 S/m (cerebrospinal fluid); the skin conductivity is set to 0.1 S/m to take into account the presence of the deep granular tissue, i.e., the dermis [39], [41]. All models were discretized with two uniform voxel meshes of resolution 1 mm and 2 mm, respectively, giving rise to the total number of voxels and unknowns, as reported in Table II. The use of different anatomical models and postures has allowed evaluating the occurrence of numerical artefacts in different exposure conditions and with different segmented models. Moreover, the use of a small anatomical model like "Roberta" (whose reduced number of voxels and unknowns corresponds to a reduced computational effort) was particularly convenient for the analysis of solver convergence reported in Section V-B.
The anatomical models are placed in close proximity to the IPT system (on the vehicle side, as sketched in Fig. 1) at a distance of 1 m from the transmitter coil axis, facing the IPT system, where the levels of the magnetic fields are maximized.
For the purpose of the accuracy assessment, we focus the attention on the exposure of a by-stander for two main reasons: it is, by now, the most common case of exposure condition analyzed in literature; moreover, it represents a good test for the analysis of the metric in presence of largely nonuniform field distribution.
The maps of the magnetic flux density magnitude (|B|) over the surfaces of the considered anatomical models are reported in Fig. 6. The corresponding maps of the raw induced E-field magnitude (|E|), also presented in Fig. 6, show how the maximum of H s is always spatially unrelated to the maximum values of the E-field. According to the ICNIRP, at the frequency of interest, the BRs (peak values) are equal to 16.2 V/m and 32.4 V/m for the general public and the workers, respectively.
For each considered formulation (T-formulation and A−φ formulation), raw electric field data are computed for 1-mm (E (1 mm) raw ) and 2-mm (E (2 mm) raw ) resolutions. The former are then averaged over 2-mm cubic volumes surrounding the voxel of interest, considering both an average independent of the adjacent tissue properties (E avg1 ) and an average involving only the voxels made of the same tissue as the reference voxel (E avg2 ). The magnitude values of the E-field raw data are also filtered using the approach for outlier removal, described in Section IV-B, obtaining the E (out) datasets.

B. Computational Efficiency and Convergence
First, the stability of the solutions has been verified by comparing the convergence of the GMRES algorithm when solving the algebraic systems that derive from the two field formulations. The analysis has been carried on with the "Roberta" model, being the computational effort lower for this model. Fig. 7 shows the behavior of the GMRES residual and the L2 norm of the approximated solution variation versus the iteration number for both formulations. The two formulations have a different convergence behavior in term of GMRES residual. Moreover, the L2 norm of the solution variation decays more rapidly for T-formulation.
The number of iterations required to reach convergence is 3976 (resp. 5050) for A−φ (resp. T) formulation, when the limit relative error of the GMRES residual is set at 10 −5 . The computational time for T formulation is higher (∼56.7 hours) with respect to the A−φ formulation (∼14.2 hours), due to the greater number of unknowns and iterations. These values refer to simulations running on a server with Intel Xeon CPU E5-2680 v2, 128 GB RAM.

C. Comparison in Terms of E-Field Values
The overall distributions of the E-field obtained with the two formulations (T and A−φ) are almost identical for all considered exposures. For this reason, Fig. 6 shows the results obtained only with the A−φ formulation for a resolution of 1 mm. However, discrepancies in terms of E-field maximum arise and depend on the considered anatomical model, exposure condition, and voxel resolution. The maximum values of the electric field magnitude, provided by the model with 1 mm resolution, are 32.7 V/m and 6.7 V/m, respectively, for the A−φ and T formulation, when the residual tolerance for GMRES is set at 10 −5 . In the next subparagraphs, the results are analyzed in detail for each considered exposure.

1) Effect of the Solver Convergence and Quantification of the Outliers:
The "Roberta" model has been adopted to analyze the effect of the GMRES solver convergence. For both formulations, Fig. 8 shows the voxel values of the magnitude of E (1 mm) raw (1 mm resolution) sorted in the ascending order. The two curves are almost coincident, apart from the last few hundreds values, where A−φ formulation shows greater outliers. The frequency distributions originated by E-field values are also similar (the diagram corresponding to the A−φ formulation is plotted in Fig. 4). The mean value of the distribution and the standard deviation are −13.09 (resp. −13.11) and 1.60 (resp. 1.64) for the A−φ (resp. T) formulation. The corresponding DP is −8.19 (resp. −8.28).  The percentage of outliers is 0.03% for both formulations. By checking the distribution of the outliers within the tissues, it results that most of them are located in the same few tissues (see Table III). In particular, the outliers that are found in skin and subcutaneous adipose tissue (SAT) are extremely thin, and in fat and bones, have low conductivities (0.043 S/m and 0.021 S/m, respectively) and presumably introduce higher discontinuities in electrical properties with respect to surrounding tissues. For this exposure, the values of the outlier data in T formulation solution are lower than in A−φ formulation solution, so, differently from what shown in Section III, the local maximum error in T formulation is presumably lower than that in A−φ formulation. Similar results are found for the other considered exposures.
A summary of the maximum values of E raw is reported in Table IV. In this table, the results obtained with the "Roberta" model exposed either to a uniform magnetic field (value equal to 20 μT, similar to the averaged IPT field) directed along the vertical axis (z-axis) or along the horizontal axis perpendicular to the thorax (y-axis) are reported. Depending on the voxel model and the distribution of the applied field (and consequently of the induced currents), the highest values of the outliers can appear with the T formulation ("Duke" standing) or the A−φ one (all other situations). This effect seems to be mainly driven by the features of the voxelized model, which may present localized areas whose voxels characteristics, interacting with the induced current paths, can enhance the appearance of artefacts. These effects can be more pronounced (see, for example, the "Roberta" model with 1 mm voxel size, either exposed to IPT field or the uniform field along the y-axis, and the "Duke" crouching model), or less evident (see, for example, the "Roberta" model exposed to the uniform field along the z-axis or the "Duke" standing model).
2) Outliers Removal: A sort of averaging effect determined by the use of larger voxel size is evident from the results presented in Table IV. As an example, for the "Roberta" model, passing from 1 to 2 mm resolution, the maximum value for the A−φ formulation reduces from 32.7 to 11 V/m (−66%), while for the T formulation reduces from 6.7 to 4.4 V/m (−33%).
Similarly, by averaging the results computed with 1 mm voxels on a 2-mm side cube (data not shown), the localized higher E-field values are only partially reduced (from 32.7 to 14 V/m for A−φ formulation and from 6.7 to 4.3 V/m for T formulation for the "Roberta" model), both when applying the averaging in each specific tissue according to [11] and when considering all tissues without distinction.
The filtering technique based on the computation of the 99th percentile of the electric field (ICNIRP recommendation) sensibly reduces the maxima (∼1 V/m, with respect to 2.4 V/m ÷ 2.5 V/m, that is a reduction of ∼60% for the "Roberta" model).
On the other hand, as already mentioned, the computation of the 99th percentile loses considerably its meaning for exposure to strongly nonuniform magnetic fields.
The raw numerical data are also filtered using the procedure discussed in Section IV-B. The results are summarized in Table IV under the column E (out) . For the "Roberta" model, the maximum value of the E-field magnitude in the 1 mm model reduces from 32.7 V/m (resp. 6.7 V/m) to 2.5 V/m (resp. 2.4 V/m) for the A−φ (resp. T) formulation. Slightly higher values are found with the 2 mm voxel size. The effect of the correction of the outliers for the raw data obtained with the A−φ and T formulations is put in evidence in Fig. 9.
In general, the application of this filtering technique allows, also for the most critical cases, a significant reduction of the higher outliers. For the "Roberta" model exposed to the IPT field, the discrepancies on the maximum detected value of the E-field (considering both voxel sizes and formulations) reduce from ∼700% to ∼37%. Similarly, for the "Duke" crouching model, the reduction is from ∼380% to ∼80% and for the "Roberta" model exposed to the uniform field along the y-axis, the reduction is from ∼660% to ∼15%.
It is worth noting that, although looking at raw data, both "Roberta" and "Duke crouching" models exceed the BRs according to the A−φ formulation, all the considered cases are sensibly below the BRs when filtered using either the 99th percentile or the more strict strategy proposed in [20].

3) Distribution of the Outliers in Different Tissues:
The comparison between the results given by the two formulations has been extended in this subparagraph to the analysis of the maximum E-field values detected within the different tissues. All results refer to the exposure of "Roberta" model to the IPT field. In Fig. 10, the maxima of the E-field peak values computed on some selected tissues is reported for the two formulations, comparing different data, and namely: the raw data computed with voxels In each diagram are reported: the raw data, computed with voxels size of 1 mm and 2 mm, the corresponding filtered values, the values obtained considering the 99th percentile evaluated tissue by tissue, and the values obtained by averaging the raw data at 1 mm on a mesh size of 2 mm. size of 1 mm and 2 mm, the corresponding filtered values, the values obtained considering the 99th percentile evaluated tissue by tissue, and the values obtained by averaging the raw data at 1 mm on a mesh size of 2 mm. The filtered values reported in the figure are obtained applying the algorithm described in Section IV-B to all voxels; the maximum values obtained on each tissue on the set of filtered data are then selected. The maximum discrepancies between the formulations arise in the skin and SAT. In particular, the T formulation seems to be more robust in terms of outliers, showing a better agreement between averaged values and raw values computed with 2 mm voxel size. Better agreement is found on internal organs. By considering the filtered data, the agreement between the maximum levels computed on each tissue by the two formulations increases and results become more stable with the voxel size.
These differences between the two formulations are almost cancelled out when saturating the maxima at the corresponding 99th percentile evaluated tissue by tissue. It is worth noting that the values given filtering the data on each single tissue are higher than the ones provided by filtering the date on the whole body, even for tissues with a large number of voxel (e.g., skin), and confirm the suspicion that 99th percentile filter could not always be a reliable metric.

VI. CONCLUSION
In this article, an analysis of the accuracy of the computational dosimetric results has been checked making reference to the human exposure to a static IPT system for the charging of electric vehicles at 85 kHz. The aim was to assess the reliability of internal field quantities (induced electric field) and, in particular, the maximum E-field value that is the quantity to be compared with the BRs for induced electric field provided by ICNIRP Guidelines 2010.
The analysis has been performed using different anatomical models (different size and posture) and exposure conditions (near-field from IPT system or uniform field) in order to investigate the main origin of local artefacts. The results confirm that the staircase interfaces, unfortunately unavoidable for voxelized human models, are responsible of most of the artefacts.
In addition, this article compares the results provided by two complementary finite-element solvers, based respectively on the use of electric scalar or vector potentials, in order to evaluate how the problem formulation could produce or enhance numerical artefacts and outliers.
The analysis has put in evidence a different behavior of the two alternative finite-element formulations for computing the induced electric field, depending on the characteristics of the voxelized anatomical human model and spatial distribution of the magnetic field.
It seems not possible to derive general conclusions in terms of the most robust method for dosimetric solutions. Nevertheless, the results obtained by the two formulations after applying the filtering techniques tend to converge.
In addition, the use of the filtering technique based on statistical considerations has demonstrated its capability to even more assess the results, removing undesirable outliers and reducing discrepancies between results obtained with different solvers to less than 10% for 1 mm anatomical model resolutions.
It must be also remarked that this filtering technique has been revisited in this article with respect to the one derived from the literature. In particular, in our implementation, some possible mismatches have been removed, demonstrating the effectiveness of this filtering technique in smoothing maximum values of induced electric fields computed in the human body model. This smoothed maximum values allow a more reliable comparison with the ICNIRP BRs in order to define as "safe" or "unsafe" a given exposure scenario related to the considered technology.
Finally, it must be pointed that the analysis has been here focused on human exposure due to IPT systems, because of the relevance of this emerging technology nowadays, but most of the considerations developed in this article could be reasonably extended to comparable exposure scenarios to low and medium frequency nonuniform magnetic fields, when reference levels cannot be adopted and compliance with the BRs is needed. Future work will be devoted to the application of the proposed methodology to the analysis of different exposure scenarios in the framework of IPT systems, including human position around and inside vehicles, human postures, and characteristics of the IPT systems (light and heavy).