In an MCR-WPTS, energy is transmitted from the transmitting side to the receiving side through electromagnetic induction [15]. The performance of the system is optimized by installing electromagnetic shielding outside the coils. Notably, the direction of the magnetic field is guided by ferrite with high-permeability, which serves to improve the coil’s mutual inductance. Meanwhile, magnetic leakage is absorbed by an aluminum plate with high conductivity. Furthermore, induced eddy currents are generated in opposite directions to reduce electromagnetic radiation.

The equivalent circuit of an MCR-WPTS equipped with two coils is presented in Fig. 2. In contrast to previous research [16], the enhanced equivalent circuit modeled in this paper accounts for the effects of the ferrite–aluminum plate composite shielding on the system. The coupling mechanism was designed as a planar structure, and the transmission medium was approximated as a homogeneous medium.

As shown in Fig. 3, the ferrite is modeled as dipole coils loaded with magnetic dipoles composed of resistance, inductance and capacitance (RLC) in series. Meanwhile, the metal aluminum plate on the receiving side, which generates the eddy currents, can be equated to a number of small inductors. It is coupled with the main coil through mutual inductance. C₁ and C₂ denote the series compensation capacitances of the transmitting and receiving sides, respectively, while C_f1 and C_f2 signify the negligible change in resonance capacitance after the addition of shielding. To maximize the transmission efficiency of the system, it is crucial to maintain the operating frequency ω as equal to the coil resonant frequency ω₀, as expressed in Eq. (1):

where L₁ and L₂ refer to the inductances of the transmitter and receiver coils, respectively, obtained through experimental measurement. Meanwhile, L_f1 and L_f2 denote changes in the self-inductance of the coils after the addition of shielding.

Based on Kirchhoff’s Voltage Law (KVL), the equation of the equivalent circuit model can be derived as follows:

Here, U_s refers to the high-frequency AC at the source, with the corresponding coil current and impedance being I and Z, respectively. M represents the mutual inductance between the transmitting and receiving coils, as formulated in Eq. (3), as noted below. Furthermore, M_Al refers to the mutual inductance between the electrically shielded aluminum plate on the receiving side and the receiving coil.

where r₁ and r₂ refer to the internal resistances of the circuit, while R_f1 and R_f2 denote the change in coil resistance after the addition of shielding, primarily indicating the hysteresis loss of the ferrite. Meanwhile, R_Al and L_Al are the resistors and inductors on the aluminum plates, respectively. M₁₂ is the mutual inductance between the original coils, and ΔM signifies the increment in mutual inductance. The equivalent load is represented as R_L.

By substituting Eq. (3) into Eq. (2) and solving the matrices, the instantaneous resonant currents on the transmitting and receiving coils were obtained. From these currents, key parameters, such as resonant voltages on the receiving coils, were derived, providing the necessary basis for the subsequent calculation of the FP on the shielding.

However, upon adding the ferrite-aluminum plate shielding, the high-frequency varying electromagnetic field induced eddy currents on the aluminum plate. This led to a significant increase in R_f and nonlinear variations in ΔM, thus complicating the numerical analysis [17].

To overcome this issue, electromagnetic simulation software COMSOL was employed to model the coil and accurately set the material parameters, allowing for the precise determination of the resonant voltage and current density of the coil.

At 85 kHz, the AC capacitive voltage to the ground and the induced voltage on the metal shielding accounted for a major part of the FP.

In Fig. 4, the working region is denoted as Ω. It is meticulously divided into Ω_I and Ω_II, where Ω_I represents the near-field region of the electromagnetic field within the shielding on both sides, and Ω_II denotes the region near the outer surface of the shielding. The boundaries of the sections, denoted as Γ, are Γ₀, Γ₁, Γ₂, and Γ₄. Γ₀ is defined as the equivalent ground boundary, Γ₁ denotes the boundary at infinity, Γ₂ indicates the receiving coil boundary, and Γ₄ is the boundary of the composite shielding layer.

The potential on the coil was uniformly distributed [18]. According to Maxwell’s equations, the relationship between the magnetic field and the field source during operation can be expressed as Eq. (4):

where n is the corresponding normal direction, ε refers to the dielectric constant of the material, and σ denotes conductivity. In an air medium, E and B are the electric field strength and magnetic induction, respectively. Furthermore, H₁, E₁, and B₁ indicate the magnetic field strength, electric field strength, and magnetic induction strength of the coil, respectively. Similarly, H₂, E₂, and B₂ are the electromagnetic parameters of electrical shielding. Lastly, the current density on the helical coil is denoted as J_s.

Upon the introduction of vector magnetic potential A and scalar function ϕ, the control equation can be rewritten as Eq. (5):

where μ refers to the material’s magnetic permeability. Assuming constant conductivity, the differential equation for the vector magnetic potential in working region Ω is given by Eq. (6):

where v is the material magnetoresistance, and J_s denotes the current density in the source region, which was equal to the current density on the coil. Notably, when the current varies sinusoidally with time while ignoring the higher harmonics of the variable, Eq. (6) can be modified into the following complex vector equation Eq. (7):

Meanwhile, the residual equation, when solved using the weighted residual method, can be formulated as Eq. (8):

where W is the vector weighting function. Furthermore, by implementing divisional integration and Green’s transform to Eq. (8), the following equation was derived as follows:

Subsequently, upon substituting and simplifying the boundary conditions, Eq. (10) was achieved:

Eq. (10) was then discretized. W represents a set of basic functions aligned with the direction of the coordinate unit vector, as represented in Eq. (11), where N represents the basis function determined by the geometry and size of the dissected cell.

By substituting the discrete Eq. (11) into the residual Eq. (8), Eq. (12) was obtained. For a simpler expression, the following variables represent the plural form:

Similarly, each of the 3n equations was obtained separately, as expressed in Eq. (13):

Eq. (14), noted below, presents these equations in matrix form:

where [S] indicates the magnetoresistive stiffness matrix, which can be defined as Eq. (15):

In this regard, the rule of construction can be expressed as Eq. (16):

where the subarray is formulated as Eqs. (17)–(19):

Here, [A] is the vector magnetic potential matrix.

Furthermore, [F] indicates the excitation matrix in Eq. (14).

Notably, the excitation of the vector magnetic potential is driven by the conduction current density. An electromagnetic analysis of the near-field region within the coupling mechanism revealed that the conduction current density was confined to the coil. Consequently, the specific vector element in the excitation matrix was obtained using Eq. (22):

In the finite element computation process, the continuous object is dissected into simple shapes to attain an approximate solution. The quality of the model is significantly influenced by the mesh, while the size of the finite element equation system is directly determined by the number of meshes. The convergence of the model is also closely related to the quality of the meshes. Additionally, higher computational accuracy implies a geometric increase in computational time. Moreover, the electrical shielding applied to the MCR-WPTS in this study is too thin compared to the whole geometry. At the same time, the use of conventional 3D profiling results in an increasing number of meshes and mesh distortions, which reduces the convergence of the numerical model. Therefore, in this study, the VHBC method was employed to optimize the [S] matrix. It is essentially a modified charge conservation method, as expressed in Eq. (23):

Upon accounting for the standard finite element matrix of Poisson’s equation, Eq. (23) can be rewritten as Eq. (24):

Notably, during the solution process, the contribution of V₂ to the global stiffness matrix is transferred entirely to the matrix elements associated with V₁ through matrix transformation. As a result, the influence of V₂ on the global stiffness matrix is represented through V₁, while the contribution of V₂ itself to [S] becomes zero after the transformation, leading to V₂ = 0. This further implies that the actual computed values of V₂ and V₁ need to be restored in the post-processing stage.

Furthermore, the Coulomb electric field near the coupling mechanism was excited by the current flowing through the coil. The boundary conditions for calculating the Coulomb potential are specified in Eq. (25):

Here, the outer boundaries of solution field Ω is denoted as Γ₀ and Γ₁. The Dirichlet boundary Γ₀ is defined as the zero potential plane, while the Neumann boundary Γ₁ is defined as the approximate infinity within the working region of the electromagnetic field. Furthermore, the potential on Γ₂ is defined as a function ϕ₀, which is numerically equal to the resonant voltage on the receiving coils. Meanwhile, Γ₄ is the electric shielding layer, the potential ϕ is the unknown quantity to be solved, FP is defined as ϕ_f, and Q_c refers to the Coulomb charge induced by the Coulomb potential on the electric shielding layer.

According to the principle of electromagnetic induction, a time-varying magnetic field induces an electric field within the conductor. In this study, the total electric field on the surface of the metal shielding was formed by a combination of this induced electric field and the Coulomb electric field generated by charge distribution. The resulting electric fields on the metal shielding were determined using Eq. (26), which provides a solution derived from the fundamental equations of electromagnetic theory.

Upon combining the determined coupling field with the boundary conditions, the results of the FP calculations were derived using Eq. (27):

where [Z] refers to the impedance matrix, m denotes the number of nodes on the Dirichlet boundary, and ϕ_i was obtained by integrating over the induced electric field. Finally, FP was calculated by adding ϕ_c and ϕ_i.

To conduct numerical calculations, a circular coupling mechanism incorporating electromagnetic shielding was designed. Single-layer coils were wound using Litz wire, maintaining perfect symmetry between the transmitting and receiving sides. Rectangular electromagnetic shields were positioned outside the coils. To satisfy the shielding requirements, the thickness of the shielding layer was constrained. The skin depth of the conductor, denoted as d, was calculated using Eq. (28):

When the metal plate’s thickness reached πd, 99.96% of the primary magnetic field strength was offset by the magnetic field, which was reversed by the eddy currents in the conductive shielding layer. Ultimately, the energy from the eddy currents was transformed into thermal energy within the conductive shielding layer.

Meanwhile, the magnetic shielding layer comprises block-shaped ferrites. Its conductivity was measured using the 4-m ratio method. Detailed parameters of the layer are provided in Table 1. The resonant voltage of the receiving coil was calculated at a frequency of 85 kHz with a supply voltage of 260 V, the results of parameters as shown in Fig. 5. During the operation of the coupling structure, electromagnetic shielding remained below the saturation magnetic flux density, ensuring constant electromagnetic characteristics.

A zero-potential condition was imposed, and the perfectly matched layer was applied as a boundary at infinity. The resonant voltage of the receiving coil, calculated as a known boundary excitation, was added and iteratively coupled using the VHBC method. The FP was calculated to be 404.75 V, as presented in Fig. 6.

During the post-processing phase, the units were standardized, and the Coulomb and induced electric fields were calculated, the results of which are illustrated in Figs. 7 and 8. The results confirm that the composite shielding structure was highly effective in suppressing the system’s magnetic field. The maximum magnetic flux density inside the shield was measured to be 0.374 T, while that at the center of the coil outside the shield was reduced significantly to 11.8 μT. Additionally, the induced electric field was successfully mitigated by the constrained magnetic field, leading to a maximum induced electric field of 140.2 V/m within the coupling structure and 11.3 V/m outside the shield. Furthermore, the electric shielding effectively limited the leakage of the Coulomb current, resulting in an external electric field of 22.31 V/m. These findings demonstrate the efficacy of the shield in minimizing both magnetic and electric field interferences, as depicted in interferences.

From Eq. (27), it is evident that FP is influenced by the parameters of the coupling mechanism. In Fig. 9, the numerical changes in FP are analyzed by varying the operating frequency of the coupling mechanism. Without accounting for frequency splitting, the FP and resonant voltage are observed to reach their peak values at 30 kHz—the FP reaches 511.6 V, while the resonant voltage peaks at 1,103.7 V. However, as the frequency increases further, these voltages gradually decrease.

Additionally, Fig. 10 illustrates the variations in FP with increasing axial and radial distances between the coils. It is observed that at an axial distance of 133 mm, the FP is measured to be 480.3 V, while the resonant voltage reaches 1,087.5 V. Meanwhile, the radial distance varies starting from zero, with the absolute value of the offset distance plotted on the horizontal axis. At a radial distance of 120 mm (absolute value), the FP rises to 675.3 V, while the resonant voltage reaches 1,365.8 V. Moreover, as the offset distance increases beyond the absolute value, both the FP and resonant voltage exhibit a rapid decline.

The above analyses indicate that the FP and resonant voltage of the receiving coil follow a similar trend—both increase initially and then undergo a decline. It is also confirmed that the structural parameters of the system influence the coupling coefficient, which subsequently affects the values of the resonant voltage and FP. Overall, these results confirm the validity of the analysis presented in Section II.

An experimental platform matching the parameters of the proposed model was constructed to validate the accuracy of the numerical calculations. The detailed parameters of the WPTS are listed in Table 2. Current sensors were installed on both the transmitting and receiving coils to accurately monitor the current flow, and a power analyzer was employed to track changes in the data. The experimental environment was set up to match the conditions used in the calculations, as shown in Fig. 11. The transmitting and receiving coils were positioned on the experimental platform. An electronic control system was employed to adjust the axial and radial distances between the coils with precision. Furthermore, an electronic load was used as the system’s load to ensure accurate testing and measurement.

In this experiment, the FP on the electrical shield was measured using a differential amplification probe. To ensure accurate and safe measurements, the differential probe and the oscilloscope ground were electrically isolated and powered independently. Notably, this isolation not only protected the oscilloscope from potential electric shocks but also extended the operational range of the measurement system [19].

The setup was designed to accurately capture the FP while minimizing interference and ensuring the safety of the measurement equipment. The FP was converted into a time-domain signal waveform, as depicted in Fig. 12, using a differential probe connected to an oscilloscope. The FP was measured to be 432 V, with the corresponding charge being −0.9 nC at a power output of 7.4 kW.

The waveform reveals the presence of harmonics, which primarily stem from electromagnetic interference issues within the public power grid. Additionally, the measured FP curve exhibits a noticeable DC component, which can be attributed to several factors. A primary contributor is the insulation resistance between the ground and the earth, which influences the selection of the absolute zero point and may introduce measurement errors. This phenomenon is illustrated in Fig. 13. In the actual experimental environment, the ground did not serve as a true zero-potential reference. Stray capacitance was present between the electrical shielding layer and the ground, while capacitance and equivalent resistance existed between the ground and the zero-potential plane. Moreover, a small grounding resistance was observed between the zero-potential plane and the point at infinity. These capacitances and resistances collectively affected the accuracy of the experimental data. It must be noted that the presence of a DC component alongside harmonics underscores the complexity of accurately measuring FP in environments where external interference and grounding issues play a significant role.

The numerical calculations of the FP were compared with the experimentally measured results, showing a discrepancy of 6.3%. This discrepancy can be explained by several factors. In addition to the errors caused by the harmonic and DC components in the FP of the WPTS, as discussed earlier, the potential contribution of the corona discharge phenomenon to the FP has not been accounted for in this study. This oversight may have further increased the observed discrepancy. The impact of corona discharge on FP will be investigated in future research to better understand and address this effect.