Open AccessArticle

Applying Finite Element Method–Dirichlet Boundary Condition Iteration to the Computation of Coupling Factors for a 3-D Human Model

Dipartimento di Ingegneria Elettrica, Elettronica e Informatica, Università di Catania, 95125 Catania, Italy

Author to whom correspondence should be addressed.

Appl. Sci. 2025, 15(2), 842; https://doi.org/10.3390/app15020842

Submission received: 6 December 2024 / Revised: 11 January 2025 / Accepted: 14 January 2025 / Published: 16 January 2025

(This article belongs to the Special Issue Electromagnetic Radiation and Human Environment)

Download

Browse Figures

Versions Notes

Abstract

The hybrid FEM-DBCI numerical method is proposed for the computation of coupling factors between time-harmonic magnetic fields and human bodies. Characteristics are highlighted which make FEM-DBCI very suitable to perform such computations. Several coil geometries are considered (circular coils) at low frequency. A simplified model of the human body is assumed as suggested by IEC standards. Details of the method are highlighted, and numerical results are also provided.

Keywords:

finite element method; integral equations; hybrid methods; eddy currents; coupling factor

1. Introduction

The increasing use of electric and electronic devices in home and work environments poses the problem of checking for the human safety of possible damages due to electromagnetic fields (EMFs), both at low and high frequencies. EMFs may cause thermal and non-thermal effects on cells, tissues, and living organisms [1]. The interactions of EMFs with the nervous system include effects on neuronal signaling, synaptic plasticity, and potential implications for neurodegenerative conditions [2]. Exposition to EMFs can lead to cardiovascular risk [3], undesired cellular and molecular effects [4], and pacemaker and cardioverter defibrillator malfunctioning [5]. Furthermore, the introduction and diffusion of new technologies, such as wireless energy transmission systems, have accentuated these problems in the case of medical implants [6].

In checking for human safety, some international organizations (IEEE, IEC, ICNIRP) have stated limits (basic restrictions) to the induced current density and electrical field in the human body [7]. However, since these limits are difficult to measure, other limits (reference levels) are stated in free space, assuming uniform fields [8]. To link basic restrictions to reference levels, several coupling factors are defined in the standards from various organizations.

Of course, there is vast scientific literature on this topic. An analytical model to establish the induced voltage in a pacemaker, considering a magnetic field in the range from 10 Hz to 30 kHz, is proposed in [9], which also reports some experimental data. In [10], measurements were performed to study safety issues due to exposures to non-uniform fields considering the coupling factor. An experimental assessment of human exposure to electromagnetic fields generated by a domestic induction cooktop was performed in [11]. The main problems related to these approaches are the low accuracy of the analytical methods, their limited application cases, and safety issues related to experimental measurements. Therefore, the use of EMF simulators is considered the best approach [12].

A 3-D human model was combined with a source model to evaluate the coupling factors in [11]. The effectiveness of the quasi-static approximation in assessing human exposure to electromagnetic fields from a wireless power transfer system has been studied in [13]. Numerical investigations to evaluate the coupling factor were reported in [14], considering the frequency range 250 Hz–10 kHz. In this study, the finite difference time domain method was used to evaluate the induced current density in an accurate human body model adopting suitable tissue parameters. The impact of electromagnetic fields on the human body is an important topic in the transportation industry, given the electrification of large vehicles and the push toward the fully a electric transportation paradigm in aircraft, ships, and vehicles [15,16]. In this context, a study on human safety issues due to EMFs in the passenger compartment of an electric vehicle was reported in [17]. The analysis focused on the internal electric field induced by the non-uniform magnetic fields generated by the cables.

In [18], different human body models were considered to evaluate the coupling factors at 50 Hz. Different coil configurations and a quasi-static electromagnetic simulation of a complex anatomical human model were considered to analyze the variation in the coupling factor versus different quantities (distance, radius, etc.). Circular coils and parallel wires with balanced currents were the two different sources of non-uniform magnetic field analyzed in [19], where numerical methods were used to assess the coupling factors under different situations in terms of non-uniform magnetic field exposure. Analyses of high-voltage power lines were led in [20], where the moment method was adopted to analyze the current induced in the body.

In this paper, we describe the use of the FEM-DBCI (Finite Element Method–Dirichlet Boundary Condition Iteration) [21,22] numerical method to perform the computation of coupling factors at low frequency by showing how this method is well suited to this aim. The most attractive property of FEM-DBCI is the possibility of leaving the various coils outside the analysis domain to be meshed, which can then be minimized. Furthermore, the same mesh can be used for all the source coils considered.

The well-known FEM-BEM (Boundary Element Method) also shares the properties of mesh minimization and invariance with respect to the change in external coils, but FEM-DBCI is simpler to implement and requires less computing time [23].

This paper is structured as follows. Section 2 briefly recalls the FEM-DBCI method. Section 3 describes how to use it for the calculation of coupling factors. Section 4 reports the results regarding some case study calculations. Finally, the authors’ conclusions follow.

2. The FEM-DBCI Method

Figure 1 depicts a typical system studied in an eddy current problem: a set of coils carry given time-harmonic source currents near some massive conductors. By assuming a quasi-static time-harmonic steady-state behavior, the problem can be mathematically formulated in terms of the following partial differential equation:

\nabla \times (μ^{- 1} \nabla \times \bar{E}) + j ω σ \bar{E} = - j ω {\bar{J}}_{s}

(1)

where

\bar{E}

is the phasor vector of the unknown electric field, μ is the magnetic permeability, σ is the electric conductivity, ω is the angular frequency, and

{\bar{J}}_{s}

is the phasor vector of the assigned source current density in the coils, assumed to be skin-effect-free.

Due to the unbounded nature of the free space around the system and in order to apply FEM, the free space is truncated by means of a closed fictitious boundary Γ_F, enclosing all the massive conductors. Optionally, some coils may be placed internally to Γ_F, but others may be left outside. On the truncation boundary Γ_F, an unknown non-homogeneous vector, the Dirichlet boundary condition, is assumed to hold:

- {\hat{n}}_{F} \times {\hat{n}}_{F} \times \bar{E} = {\bar{E}}_{F}

(2)

where

{\hat{n}}_{F}

is the outward unit vector normal to Γ_F and

{\bar{E}}_{F}

is the unknown component of the electric field tangent to the fictitious boundary.

The bounded domain inside Γ_F is discretized by means of tetrahedral edge elements [21]. In each tetrahedron, the electric field is developed as

\bar{E} = \sum_{m = 1}^{6} E_{m} {\bar{α}}_{m}

(3)

where

{\bar{α}}_{m}

are the non-dimensional vector shape functions:

{\bar{α}}_{m} = L_{m} (ξ_{1 m} \nabla ξ_{2 m} - ξ_{2 m} \nabla ξ_{1 m})

(4)

with ξ_1m and ξ_2m local coordinates (in the tetrahedron) associated with the first and second node, respectively, of the m-th edge e_m, whose length is L_m; the scalar quantity E_m (in V/m) is the m-th expansion coefficient of the electric field and is given by the mean value of the tangent component of the field along the edge:

E_{m} = \frac{1}{L_{m}} \int_{e_{m}} \bar{E} \cdot {\hat{t}}_{m} d s

(5)

where

{\hat{t}}_{m}

is the unit vector along the m-th edge.

By using the vector shape functions

{\bar{α}}_{m}

as weighting functions, the Galerkin method leads to the following matrix equation:

A e = b_{0} - A_{F} e_{F}

(6)

where A and A_F are sparse matrices depending on geometry and materials, e and e_F are the unknown vectors of the field expansion coefficients E_m for the internal and boundary edges, respectively, and b₀ is a vector due to the internal source currents, if there are any.

Another equation relating e to e_F can be derived by means of the following representation formula:

\frac{γ}{4 π} \bar{E} (P_{F}) = {\bar{E}}_{e x t} (P_{F}) + \frac{1}{4 π} ∯_{Γ} (\frac{1}{r} \hat{n} \times \nabla \times \bar{E} + (\hat{n} \times \bar{E}) \times \nabla \frac{1}{r} + \hat{n} \cdot \bar{E} \nabla \frac{1}{r}) d S

(7)

where P_F is a point lying on Γ_F, Γ is another closed surface that includes the system (see Figure 1), r is the distance between P_F and the generic point P on Γ,

{\bar{E}}_{e x t}

is the electric field due to the coil currents external to Γ_F, and γ is a solid angle that is equal to 4π if P_F is external to Γ; if P_F lies on Γ, γ is the solid angle at P_F of the unbounded free space external to Γ. The physical meaning of Equation (7) is that if a vector field and its rotation are known on a closed surface Γ, then the field is completely computable in the entire space outside that surface.

The closed surface Γ is conveniently selected as constituted of the triangular faces of the tetrahedral elements. If the surface Γ is selected as strictly enclosed by Γ_F, the integrand function in Equation (7) is regular and the method is referred to as FEM-DBCI [21]. Otherwise, if the surface Γ is selected as coincident with Γ_F, the integrand function in (7) is singular and the method is called FEM-SDBCI (Singular DBCI) [22].

By substituting (7) in (5), the following integral equation is obtained:

\begin{array}{l} \frac{γ}{4 π} E_{m} = \frac{γ}{4 π L_{m}} \int_{e_{m}} \bar{E} (P_{F}) \cdot {\hat{t}}_{m} d s = \frac{1}{L_{m}} \int_{e_{m}} {\bar{E}}_{e x t} (P_{F}) \cdot {\hat{t}}_{m} d s + \\ + \frac{1}{4 π L_{m}} \sum_{k} \iint_{T_{k}} \int_{e_{m}} (\frac{1}{r} \hat{n} \times \nabla \times \bar{E} + (\hat{n} \times \bar{E}) \times \nabla \frac{1}{r} + \hat{n} \cdot \bar{E} \nabla \frac{1}{r}) \cdot {\hat{t}}_{m} d s d S \end{array}

(8)

where T_k is the k-th triangular patch of the integration surface Γ, resulting from the tetrahedral meshing of the bounded domain. In each patch, the electrical field is developed as in (3), where the six edges are those of the corresponding tetrahedron to which the triangular patch belongs.

Finally, by developing the electrical field as in (3) and by applying (8) to all the edges e_m of the fictitious boundary, we have the following:

H e_{F} = e_{e x t} + G e

(9)

where e_ext is the vector of the mean tangent components of the electrical field along the edges of the fictitious boundary Γ_F due to the external coils, and H and G are dense matrices collecting the geometrical coefficients in (8) relative to the edges on Γ_F and to the internal edges, respectively. Note that not all internal edges appear in (8), but only those of the tetrahedra external to Γ and with a face lying on it. So, the matrix H is square by construction, while G, in general, is rectangular. Note that in the case of FEM-DBCI, matrix H simplifies to the identity matrix and Equation (9) becomes explicit with respect to vector e_F.

Gauss quadratures are used to evaluate both the double integral on the triangle T_k and the line integral on the edge e_m in (8). In order to obtain a good tradeoff between accuracy and computing speed, the following rule is recommended. If L_max is the length of the longest edge of the triangle T_k on Γ, L_m is the length of the edge on the fictitious boundary, L = max (L_max, L_m), and d is the distance between their centers, then for L/d ≤ 0.2, a one-point quadrature is used on both the triangle and the edge; for 0.2 < L/d ≤ 1.1, three Gauss points are used on the triangle and two points on the edge; otherwise, six points are used on the triangle and three points on the edge. In the case of FEM-SDBCI, the singularities arising in the integrand function in (8) are overcome by means of analytical formulas [24]. Some details are given in [22].

Combining the matrix Equations (6) and of (9), the global linear system of the FEM-DBCI (or FEM-SDBCI) method is obtained:

[\begin{matrix} A & A_{F} \\ - G & H \end{matrix}] [\begin{matrix} e \\ e_{F} \end{matrix}] = [\begin{matrix} b_{0} \\ e_{e x t} \end{matrix}]

(10)

The system matrix in (10) is partly sparse, partly full, and non-symmetric. Unfortunately, efficient solvers do not exist for this kind of linear algebraic system. However, a good solving strategy should treat the sparse subsystem coming from the FEM equations and the dense subsystem coming from the integral equations in different ways. In particular, this strategy should reduce the number of multiplications of the dense submatrices by vectors to a minimum and should fully exploit the sparsity and symmetry of matrix A. A strategy that satisfies the above requirements is the following two-block Gauss–Seidel algorithm [22]:

(1): select arbitrarily a first guess for e_F, for example, e_F = 0;
(2): solve (6) for e by means of the Complex Conjugate Gradient (CoCG) [25];
(3): in the case of FEM-SDBCI decompose the square matrix H is into the matrices L and U; then Equation (9) is solved for e_F; in the case of FEM-DBCI the vector e_F is computed directly as $e_{F} = e_{e x t} + G e$ ;
(4): to measure the distance between the new guess for e_F and the old one, the convergence indicator η is computed:

$η = 100 \frac{{‖e_{F}^{n e w} - e_{F}^{o l d}‖}_{2}}{{‖e_{F}^{n e w}‖}_{2}}$

(11)

convergence is assumed if η < η₀, where η₀ is a user-selected end-iteration tolerance;
(5): if convergence is not reached, go to step 2 by assuming a relaxed new guess for e_F such as:

$e_{F}^{(n)} = λ e_{F}^{n e w} + (1 - λ) e_{F}^{o l d}$

(12)

where λ is the relaxation coefficient.

This iterative algorithm exhibits the following characteristics:

(a): Since in each iteration, the first guess for the CoCG is the solution obtained in the previous iteration, the various solutions of the FEMl Equation (6) get faster as the iteration proceeds.
(b): In the case of FEM-SDBCI, the LU decomposition of matrix H is performed only once at the beginning of the iterative procedure; round-off errors in the LU decomposition are treated by selecting an appropriate (double precision) accuracy in computing and storing matrices H and G.
(c): The whole iterative procedure is convergent if an appropriate relaxation coefficient λ is selected, which, unfortunately, is not known a priori; divergence may occur if an unsuitable coefficient λ is used.
(d): Consequently, the integral Equation (9) is used only a few times, compared to its use in an iterative CG-like solver for the whole non-symmetric system (10).

Another solving strategy has been devised by considering the following reduced algebraic system:

M e_{F} = b

(13)

where

M = H + G A^{- 1} A_{F}

(14)

b = e_{e x t} + G A^{- 1} b_{0}

(15)

Of course, the definitions of matrix M in (14) and vector b in (15) are only formal because the inverse matrix A⁻¹ cannot be really obtained when the number of unknowns is large. However, these arrays can be used, since vector b is directly computable, while matrix M can be virtually used to perform matrix multiplications by vectors, as explained below.

In order to compute vector b, the following must be performed: solve the FEM Equation (6) with e_F = 0 by means of the CoCG solver to obtain e = A⁻¹b₀, and then compute b = e_ext + Ge. Similarly, in order to compute the product Me_F of matrix M by a given vector e_F, solve FEM Equation (6) with b₀ = 0 using the CoCG solver to obtain e = A⁻¹A_Fe_F, and then compute Me_F = He_F + Ge.

Note that the iterative algorithm described above can be seen as the application of the Richardson method to solve (13). However, since this solution method is known to be weak and possibly non-convergent, one should look for better solvers.

Due to the virtual availability of matrix M, it is possible to use various non-stationary solvers, and in particular, solvers that are polynomial accelerations of the Richardson method for non-symmetric matrices. Among these, the Biconjugate Gradient (BiCG), the Quasi Minimal Residual (QMR), the Conjugate Gradient Squared (CGS), the Biconjugate Gradient Stabilized (BiCGstab), and the Generalized Minimal Residual (GMRES) are the most popular [25].

Since, concerning system (13), the matrix–vector multiplications are much more expensive than in a system where the coefficient matrix is directly available, it is crucial to reduce the number of iterations. Then, GMRES should be preferred, because it minimizes the number of iterations by performing a true minimization of the residual [26]. Moreover, note that the residual can be computed directly with the approximate solution (but this requires a further matrix–vector multiplication) or by using the orthonormal basis of the Krylov subspace [26]. The latter option is preferable. In fact, in general, the greater drawbacks of GMRES are the computing time and memory required to compute and store the orthonormal basis, which increases linearly with the number of iterations, but in this context, the unknowns in system (13) are the values of the electric field along the edges of the fictitious boundary, the number of which is relatively low compared to the total number of edges.

Another important point regards the GMRES restating procedure. It is convenient to use long restarts, which generally result in a full GMRES due to the quick convergence characteristic of the simple iterative procedure. The fact that the relaxed iterative procedure converges with a suitable choice of a positive relaxation parameter λ indicates that the eigenvalues of matrix M have positive real parts, and this assures that GMRES converges to the true solution even with a very short restarting parameter.

Comparing the GMRES solution algorithm with the simple iterative one, it is noteworthy that the GMRES solution does not require the LU decomposition of the H matrix. This is a great advantage, especially in the case of problems with a large number of unknowns. Conversely, the various solutions of the FEM equations by means of the CoCG-solver are not related to each other, so the number of CoCG steps does not decrease as the iteration proceeds.

3. Applying FEM-DBCI to Computation of Coupling Factors

The FEM-DBCI and FEM-SDBCI methods described in the previous section appear very suitable for the computation of coupling factors, as can be understood by considering the following points.

(a): The finite element mesh is composed of tetrahedra which are well suited to describe complex bodies, even the most realistic and complex human models, with several types of tissues.
(b): The tetrahedral mesh can be minimized, since there is no need to discretize any free space in the FEM-SDBCI by selecting the integration surface Γ as the skin of the human body and the fictitious boundary Γ_F to coincide with it. In the case of FEM-DBCI, it is sufficient to mesh only a thin layer of air by selecting Γ as before and Γ_F at a small distance from the human skin.
(c): These reduced meshes allow us to place the source coils outside the analysis domain, thus vector b₀ in (6) vanishes and the various matrices A, A_F, H, and G do not change from one numerical analysis to another by modifying the external coil. Only the e_ext vector depends on the actual coil type and position; the entries of this array are computed by means of

$\bar{E} (P_{F}) = - j ω \bar{A} (P_{F})$

(16)

where $\bar{A} (P_{F})$ is the magnetic vector potential, which can be computed analytically if one considers simple coil geometries.
(d): Note that the A and A_F matrices in the FEM equation are complex, while H and G in the integral equation are real; in fact, these matrices depend only on the geometry of the surfaces Γ_F and Γ and on the mesh; hence the matrices H and G do not change with the (low) frequency of the source currents. Due to the low conductivity of the human model, the skin depth is very large, so it is possible to use a tetrahedral finite element mesh with an average edge size of about 3 cm, without modifications up to frequencies of 10 megahertz.
(e): In post-processing, the magnitudes of the electrical field E, and hence of the current density J = σE, are computed at the barycenters of the various tetrahedra in the mesh and also at the barycenters of the various triangular faces on the skin of the human model:

$E = |\bar{E}| = |\sum_{m = 1}^{6} E_{m} {\bar{α}}_{m} (P_{n})|$

(17)

where P_n denotes the generic barycenter. Note that the heating power density (in W/m³) is obtained as ½σE².
(f): The computation of the flux density B, if needed, can be performed in post-processing on a per-element basis:

\bar{B} = - \frac{1}{j ω} c u r l \bar{E} = - \frac{1}{j ω} \sum_{m = 1}^{6} E_{m} c u r l {\bar{α}}_{m}

(18)

Since the vector shape functions α_m vary linearly in space, the flux density in (18) is constant in each tetrahedron. As will be explained in the next section, the magnitude of the magnetic field B₀ must be computed in the tetrahdron of the human model nearest to the coil.

(g): For the solution of the hybrid global systems, the use of the GMRES, as described in the previous section, is more suitable since its convergence characteristics appear more robust than that of the simple iteration. The restating parameter has been set to 10.

4. Numerical Computations

A simple 3-D human model (head and torso) is used for the computations as suggested in the IEC 62311 document [8]. The electrical conductivity σ is set to 0.2 S/m, while free space permeability μ₀ is assumed. The lowest point of the model is located at the origin, while the z-axis is its rotational symmetry axis. The height of the human model is set to h₀ = 152.8 cm. An ideal circular coil is considered with its axis parallel to the z-axis, radius R_c = 16 cm, placed at z_c = 100 cm from the soil (the xy plane). Its distance d from the model varies from 10 to 30 cm. A time-harmonic source current is assumed to flow in the coil, with intensity I = 1 A and frequency f = 50 Hz. Due to symmetry reasons, only half the human model is meshed by assuming a homogeneous Dirichlet boundary condition on the xz symmetry plane. Figure 2 shows the adopted tetrahedral mesh (3067 nodes, 12,882 tetrahedra, and 17,764 edges) as seen frontally from a viewpoint in the negative y-axis. The average length of the edges is about 3 cm.

Indicating with C ≡ (x_C, 0, z_C) the center of the circular coil, the three Cartesian components of the magnetic vector potential at a point P = (x, y, z) in free space are given by

A_{x} (P) = - A_{φ} \sin (φ), A_{y} (P) = A_{φ} \cos (φ), A_{z} (P) = 0

(19)

where

A_{φ} = \frac{μ_{0} I}{π k} \sqrt{\frac{R}{ρ}} [(1 - k^{2} / 2) E (k) - F (k)]

(20)

ρ = \sqrt{{(x - x_{C})}^{2} + y^{2}}

(21)

k = \sqrt{\frac{4 R ρ}{{(R + ρ)}^{2} + {(z - z_{C})}^{2}}}

(22)

and E(k) and F(k) are the complete elliptic integrals of the first and second kinds, respectively, of modulus k.

In order to compute the generic entry E_m of the vector e_F, the Gauss quadrature is used for integration along the edge e_m lying on Γ_F:

E_{m} = \frac{1}{L_{m}} \int_{e_{m}} \bar{E} \cdot {\hat{t}}_{m} d s = - \frac{j ω}{L_{m}} \sum_{i = 1}^{N g} w_{i} (A_{x} (P_{i}) e_{x} + A_{y} (P_{i}) e_{y})

(23)

where N_g is the number of Gauss points P_i (i = 1,…,N_g) on the edge e_m, w_i (i = 1,…,N_g) are their weights, and e_x and e_y are the x- and y-components of the edge. The number N_g is selected as follows: if d_m is the distance of the barycenter B_m of edge e_m from the coil point P₀≡(x_C − R_C,0,z_C) nearest to the model, then for L_m/d ≤ 0.2 a one-point quadrature is used; for 0.2 < L_m/d ≤ 1.1, two Gauss points are used; otherwise, three points are used on the edge.

The FEM-SDBCI method was used to analyze this system, so that the human skin is at the same time both the fictitious boundary Γ_F and the integration surface Γ, and the various coils are left outside the domain of the FEM analysis.

The solution has been obtained by means of the GMRES solver as described in Section 2; only seven iterations were performed to reach convergence with an end-iteration tolerance of η₀ = 0.01 percent. From the solution, the modulus of the electrical field was evaluated at each barycenter of the tetrahedra inside the body and of the triangles on the skin, and the maximum value E_max is obtained. From this, we obtain the maximum value J_max = σE_max of the current density in order to compute the coupling factor.

K = \frac{J_{\max}}{J_{\max}^{(u)}}

(24)

Also, the maximum value

J_{\max}^{(u)}

of the current density needs to be computed. This value must be evaluated in the situation in which the human body is immersed in a uniform external magnetic field whose intensity is equal to that of previous analysis at the point at which the human skin is nearest to the coil.

To compute

J_{\max}^{(u)}

, another numerical analysis must be conducted. In order to obtain a quasi uniform magnetic field, the source may be constituted by a couple of horizontal circular coils, with a large radius R₀, centered at the z-axis and located at heights (h₀ + R₀)/2 and (h₀ − R₀)/2 from the xy plane, where h₀ is the eight of the human model. The two coils carry the same sinusoidal current of I = 1 A. Of course, the use of the FEM-DBCI or FEM-SDBCI methods allows us to utilize the same mesh and hence the same global algebraic system of the previous analysis, except for the known term of the integral Equation (9). Once this new solution is obtained, one can compute the maximum value of the current density

J_{\max}^{(1)}

in the human model and the flux density

B_{0}^{(1)}

in the same tetrahedron where the flux density B₀ was previously evaluated. Then, due to the linearity of the material of the human model,

J_{\max}^{(u)}

in (24) is evaluated as

J_{\max}^{(u)} = \frac{B_{0}}{B_{0}^{(1)}} J_{\max}^{(1)}

(25)

In Figure 3, the behavior of the coupling factor K is given versus the distance d. As can be seen, the values obtained are in good agreement with those reported in [18].

In Figure 4, the heating power density is drawn on the skin of the human model in front of the circular coil, for two different placements of the circular coil at distances of 10 cm and 20 cm.

5. Conclusions

The hybrid numerical methods FEM-DBCI and its variant FEM-SDBCI were proposed to evaluate the coupling factors between the human body and the external coils in the case of low frequency time-harmonic currents. The more attractive feature of these methods is the possibility of leaving the various coils outside the FEM mesh of the domain of analysis.

In this way, the same finite element mesh can be utilized for all the positions and kinds of the source coils to be considered. Moreover, the free space mesh around the human model is reduced to a minimum in FEM-DBCI, or really to zero in FEM-SDBCI. Note that several commercial codes do not have effective treatments for open boundary domains, so much of the finite element mesh is in free space, since the truncation boundary (where a homogeneous boundary condition is very often imposed) must be placed far from the system core in order to obtain acceptable results. In addition, a further difficulty arises from the need to mesh the source coil with finite elements of very different sizes than those of the human body.

Further research plans to use more complicated human models and to develop an edge version of the non-standard FEM-BEM method proposed in [27] for nodal finite elements.

The FEM-DBCI and FEM-SDBCI methods were implemented in ELFIN, an FEM code developed by the authors for research in computational electromagnetics [28].

Author Contributions

Conceptualization, G.A., S.A. and N.S.; Software, S.A., S.A.R. and N.S. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data is contained within the article.

Conflicts of Interest

The authors declare no conflict of interest.

References

Wust, P.; Kortüm, B.; Strauss, U.; Nadobny, J.; Zschaeck, S.; Beck, M.; Stein, U.; Ghadjar, P. Non-thermal effects of radiofrequency electromagnetic fields. Sci. Rep. 2020, 10, 13488. [Google Scholar] [CrossRef]
van der Meer, J.N.; Eisma, Y.B.; Meester, R.; Jacobs, M.; Nederveen, A.J. Effects of mobile phone electromagnetic fields on brain waves in healthy volunteers. Sci. Rep. 2023, 13, 21758. [Google Scholar] [CrossRef]
Vangelova, K.; Deyanov, C.; Israel, M. Cardiovascular risk in operators under radiofrequency electromagnetic radiation. Int. J. Hyg. Environ. Health 2006, 209, 133–138. [Google Scholar] [CrossRef] [PubMed]
Lai, H.; Levitt, B.B. Cellular and molecular effects of non-ionizing electromagnetic fields. Rev. Environ. Health 2023, 39, 519–529. [Google Scholar] [CrossRef]
Upadhyay, S.; Upadhya, A.; Salehi, W.; Gupta, G. The medical aspects of EMI effect on patients implanted with pacemakers. Mater. Today Proc. 2021, 45, 5243–5248. [Google Scholar] [CrossRef]
Christ, A.; Douglas, M.; Nadakuduti, J.; Kuster, N. Assessing human exposure to electromagnetic fields from wireless power transmission systems. Proc. IEEE 2013, 101, 1482–1493. [Google Scholar] [CrossRef]
C95.1-2019 IEEE; Standard for Safety Levels with Respect to Human Exposure to Electric, Magnetic, and Electromagnetic Fields, 0 Hz to 300 GHz. (Revision of IEEE Std C95.1-2005/ Incorporates IEEE Std C95.1-2019/Cor 1-2019). IEEE: New York, NY, USA, 2019. [CrossRef]
I.S. EN IEC 62311:2020; Assessment of Electronic and Electrical Equipment Related to Human Exposure Restrictions for Electromagnetic Fields (0 Hz-300 GHz). IEC: Geneva, Switzerland, 2007.
Hille, S.; Eichhorn, K.F.; Gonschorek, K.-H. Determination of the interference voltage in implantable medical devices with bipolar electrodes. In Proceedings of the International Symposium on Electromagnetic Compatibility—EMC Europe, Hamburg, Germany, 8–12 September 2008. [Google Scholar] [CrossRef]
Sunohara, T.; Hirata, A.; Laakso, I.; de Santis, V.; Onishi, T. Evaluation of nonuniform field exposures with coupling factors. Phys. Med. Biol. 2015, 60, 8129–8140. [Google Scholar] [CrossRef]
Jung, K.-J.; Byun, J.-K. The human exposure assessment of magnetic field from an induction cooktop using coupling factor based on measurement data. J. Magn. 2018, 23, 473–479. [Google Scholar] [CrossRef]
Ahn, J.; Hong, S.-E.; Kim, H.; Song, K.; Choi, H.-D.; Ahn, S. Improved calculation method of coupling factors for low-frequency wireless power transfer systems. Int. J. Environ. Res. Public Health 2021, 19, 44. [Google Scholar] [CrossRef]
Laakso, I.; Shimamoto, T.; Hirata, A.; Feliziani, M. Quasistatic approximation for exposure assessment of wireless power transfer. IEICE Trans. Commun. 2015, E98B, 1156–1163. [Google Scholar] [CrossRef]
Gustrau, F.; Bahr, A.; Rittweger, M.; Goltz, S.; Eggert, S. Simulation of induced current densities in the human body at industrial induction heating frequencies. EEE Trans. Electromagn. Compat. 1999, 41 Pt 2, 480–486. [Google Scholar] [CrossRef]
Barzkar, A.; Ghassemi, M. Components of Electrical Power Systems in More and All-Electric Aircraft: A Review. IEEE Trans. Transp. Electrification 2022, 8, 4037–4053. [Google Scholar] [CrossRef]
Alonso-Cepeda, A.; Villena-Ruiz, R.; Honrubia-Escribano, A.; Gómez-Lázaro, E. A veview on electric vehicles for holistic robust integration in cities: History, legislation, meta-analysis of technology and grid impact. Appl. Sci. 2024, 14, 7147. [Google Scholar] [CrossRef]
Hakuta, Y.; Watanabe, T.; Takenaka, T.; Ito, T.; Hirata, A. Safety standard compliance of human exposure from vehicle cables using coupling factors in the frequency range of 0.3–400 kHz. IEEE Trans. Electromagn. Compat. 2020, 63, 313–318. [Google Scholar] [CrossRef]
Jung, K.-J.; Shim, J.-H.; Choi, M.-S.; Byun, J.-K. Non-uniform magnetic field exposure assessment using coupling factors based on 3-D anatomical human model. EEE Trans. Magn. 2017, 54, 5000104. [Google Scholar] [CrossRef]
Shim, J.-H.; Choi, M.-S.; Jung, K.-J.; Kwon, J.-H.; Byun, J.-K. Comparative study of coupling factors for assessment of low-frequency magnetic field exposure. J. Magn. 2016, 21, 516–523. [Google Scholar] [CrossRef]
Chen, H.-Y.; Peng, C.-K.; Lin, D.-P. Currents induced in the human body by ELF electric fields inside houses located under power lines. Radio Sci. 1999, 34, 1013–1024. [Google Scholar] [CrossRef]
Aiello, G.; Alfonzetti, S.; Dilettoso, E. Finite element solution of eddy current problems in unbounded domains by means of the hybrid FEM−DBCI method. IEEE Trans. Magn. 2003, 39, 1409–1412. [Google Scholar] [CrossRef]
Aiello, G.; Alfonzetti, S.; Dilettoso, E.; Salerno, N. Eddy current computation by the FEM-SDBCI method. IEEE Trans. Magn. 2016, 52, 7400804. [Google Scholar] [CrossRef]
Aiello, G.; Alfonzetti, S.; Borzì, G.; Rizzo, S.A.; Salerno, N. A comparison between hybrid methods: FEM-BEM versus FEM-DBCI. COMPEL 2013, 32, 1901–1911. [Google Scholar] [CrossRef]
Graglia, R.D. On the numerical integration of the linear shape functions times the 3-D Green’s function or its gradient on a plane triangle. IEEE Trans. Antennas Propag. 1993, 41, 1448–1455. [Google Scholar] [CrossRef]
Golub, G.H.; Van Loan, C.F. Matrix Computations; The John Hopkins University Press: London, UK, 1996. [Google Scholar]
Saad, Y.; Schultz, M.H. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J. Sci. Stat. Comput. 1986, 7, 856–869. [Google Scholar] [CrossRef]
Alfonzetti, S.; Salerno, N. A non-standard family of boundary elements for the hybrid FEM-BEM method. IEEE Trans. Magn. 2009, 45, 1312–1315. [Google Scholar] [CrossRef]
Aiello, G.; Alfonzetti, S.; Borzì, G.; Salerno, N. An Overview of the ELFIN Code for Finite Element Research in Electrical Engineering. In Software for Electrical Engineering: Analysis and Design IV; Konrad, A., Brebbia, C., Eds.; WIT Press: Southampton, UK, 1999. [Google Scholar]

Figure 1. Massive eddy current conductors (characterized by conductivity σ and permeability μ) and distributed internal source coils (characterized by the current density J) are enclosed by the fictitious truncation boundary Γ_F. Other lumped and distributed coils are left outside (external coils).

Figure 2. Tetrahedral mesh of the human model and its enlargement in the head part.

Figure 3. Behavior of the coupling factor K versus the distance d.

Figure 4. Distribution of the heating power density on the skin of the human model at two different distances (d = 10 cm and d = 20 cm) of the coil from the skin of the human model.

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2025 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Aiello, G.; Alfonzetti, S.; Rizzo, S.A.; Salerno, N. Applying Finite Element Method–Dirichlet Boundary Condition Iteration to the Computation of Coupling Factors for a 3-D Human Model. Appl. Sci. 2025, 15, 842. https://doi.org/10.3390/app15020842

AMA Style

Aiello G, Alfonzetti S, Rizzo SA, Salerno N. Applying Finite Element Method–Dirichlet Boundary Condition Iteration to the Computation of Coupling Factors for a 3-D Human Model. Applied Sciences. 2025; 15(2):842. https://doi.org/10.3390/app15020842

Chicago/Turabian Style

Aiello, Giovanni, Salvatore Alfonzetti, Santi Agatino Rizzo, and Nunzio Salerno. 2025. "Applying Finite Element Method–Dirichlet Boundary Condition Iteration to the Computation of Coupling Factors for a 3-D Human Model" Applied Sciences 15, no. 2: 842. https://doi.org/10.3390/app15020842

APA Style

Aiello, G., Alfonzetti, S., Rizzo, S. A., & Salerno, N. (2025). Applying Finite Element Method–Dirichlet Boundary Condition Iteration to the Computation of Coupling Factors for a 3-D Human Model. Applied Sciences, 15(2), 842. https://doi.org/10.3390/app15020842

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Applying Finite Element Method–Dirichlet Boundary Condition Iteration to the Computation of Coupling Factors for a 3-D Human Model

Abstract

1. Introduction

2. The FEM-DBCI Method

3. Applying FEM-DBCI to Computation of Coupling Factors

4. Numerical Computations

5. Conclusions

Author Contributions

Funding

Institutional Review Board Statement

Informed Consent Statement

Data Availability Statement

Conflicts of Interest

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI