Design Optimization of Inductive Power Transfer Systems Considering Bifurcation and Equivalent AC Resistance for Spiral Coils

This paper presents a design optimization algorithm for series-series compensated Inductive Power Transfer (IPT) system based on flat spiral coils, considering bifurcation phenomenon and AC equivalent resistance of the coils. Moreover, it finds the best values in the areas where the IPT system operates in Zero Voltage Switching (ZVS) condition. Equivalent AC resistance of spiral coils is modeled based on eddy currents simulations using Finite Element Method (FEM) and Maxwell simulator. Based on the FEM simulations, a new approximation method using separation of variables is proposed as a function of spiral coil’s main parameters. This paper shows that this model accurately derive equivalent resistance of the coils for a specific strand diameter, with almost 95% accuracy. Based on the proposed algorithm, several IPT systems are optimized in MATLAB software using Genetic Algorithm (GA). Finally, the proposed method has been verified using a laboratory prototype with output power of about 200 watts and operating frequency of about 85 kHz.


INDEX TERMS
Inductive Power Transfer (IPT) is an efficient and reliable solution for charging of Electric Vehicles (EVs) and mobile devices, which removes heavy and troublesome cables. Moreover, IPT systems have unique advantages such as being not sensitive to environmental conditions, better maintenance and highly reliable [1]- [5]. Due to these advantages, IPT systems have wide applications for different power levels.
High-power applications of IPT systems include wireless charging of EVs and public transportation systems [6]- [8]. Low-power applications of IPT system include wireless charging of biomedical implants [9], [10] and wireless charger for mobiles devices [2], [11]. Inductive chargers are based on magnetic coupling driven by a high frequency alternating source. The power sources are usually based on resonant inverters, which are classified according to their resonant tank topologies. Common types of topologies include Series-Series (SS), Series-Parallel (SP), Parallel-Series (PS) and Parallel-Parallel (PP) topologies. Several articles analyze and compare the behavior of the topologies [12]- [14]. Among them, the SS topology has a desired characteristic that makes the compensating capacitors independent of the load and magnetic coupling coefficient variations.
Despite of wide applications for IPT systems, there are challenges associated with their design optimization. One of the main concerns is how to design the coils for a desired power with maximum efficiency considering the all physical parameters, i.e. airgap, number of strands, number of turns, etc. Moreover, the design procedure gets more complicated considering Zero Voltage Switching (ZVS), bifurcation phenomenon and AC equivalent resistance of the coils. Previous researches have no general optimization and some of them present a design procedure even without considering AC resistance [15], [16].
AC resistance of a coil is significantly larger than DC resistance under higher operating frequencies [17]- [20]. Furthermore, previous works have no discussion about ZVS operating area and avoid bifurcation phenomenon in their designs [20]- [23], which have less generality.
In the proposed algorithm, parameters of the IPT system including values of the primary and secondary inductances, number of strands, space between turns and dimensions of the spiral coils are considered as degrees of freedom. Moreover, achieving ZVS area, bifurcation phenomenon and modeling of AC resistance of Litz wires are considered in the proposed algorithm.
In order to obtain AC equivalent resistances, eddy current analysis are done using Finite Element Method (FEM). Based on the FEM simulations, new polynomial equations for modeling of the AC resistances are devised as a function of the spiral coil parameters, using a simple separated variables method, with almost 95% accuracy. The results show that the proposed method can be executed for each Standard Wire Gauge (SWG).
In this paper, the proposed modeling is done for operating frequency of about 85 kHz, following SAE J2954 recommendations and some recent works [4], [24], [25]. This paper shows that the proposed algorithm based on the AC resistance model is a universal solution for design optimization of SS compensated IPT systems. This paper is organized as follows: Section II presents modeling of the IPT system, flat spiral coils and the proposed AC resistance modeling. Section III presents the optimization algorithm and its results. Section IV and V present the experimental results of the laboratory prototype and the main conclusion of the paper, respectively. Fig.1 shows the main blocks of the IPT system, such as transmitter and receiver coils, compensating capacitors, and inverter and rectifier units. The IPT system can be modeled based on First Harmonic Approximation (FHA), which is shown in Fig. 2. The DC source and half-bridge inverter are modeled by an alternative voltage source, V p , at the fundamental frequency. It should be mentioned that the proposed design optimization can also be applied for full bridge inverters. In addition, the primary and secondary capacitances, C p and C s , primary and secondary selfinductances, L p and L s , mutual inductance, M , and equivalent load resistance, R e , are shown based on a T-model of the coupled inductors.  Assuming the secondary voltage, V s , is the reference phasor, the primary voltage and current, V p and I p , are derived by (1) and (2), respectively.

A. MODELING OF THE IPT SYSTEM
where, θ and φ are the input impedance phase angle and phase difference between the primary and secondary voltages, respectively. First harmonic equations of the secondary voltage and current, V s and I s , are derived as follows, which are considered as the reference phasor: At the beginning, the Constant Current (CC) mode is utilized and battery voltage increases up to its nominal voltage [3], [5], and [26]. Hence, the maximum, power of the charger occurs at the end of CC mode and nominal equivalent load resistance seen from secondary coil, R e , is obtained using (3) and (4), as follows: where, R L is obtained from the nominal battery voltage, V Bat , and current, I Bat , at the end of CC mode. The natural angular frequency of the IPT system is derived from (7): By the phasor analysis of the circuit presented in Fig. 2, the relationship between the primary and secondary currents is specified in (8): Total impedance seen from the primary, Z in and its phase, θ, is obtained from (9), as a function of angular switching frequency, ω s = 2π f s .
where, θ is the input impedance phase angle. According to equations (1) and (9), the amplitude of the primary current is derived as follows: Regarding (6) up to (9), the input power, P in , output power, P out , and efficiency, η, are obtained by: The output power and efficiency equations are the key equations to optimize the IPT system.

B. MODELING OF FLAT SPIRAL COILS
The spiral geometry is more common and proper for static IPT systems [2], [7] and [12]. Hence, this paper uses this geometry for the transmitter and receiver coils. Modeling of spiral coils is important for designing. According to the formula presented by Wheeler [27], to calculate the outer diameter, the following equation can be used: where, D i is the inner diameter, S is the space between adjacent turns, D w is the diameter of Litz wire and N is the turn number of coil. Regarding [27], equation (15) can be used to calculate the inductance in µH, while the dimensions are in inch.
The coupling between primary and secondary coils is obtained from (16): According to Ohm's law, the DC resistance, R dc , of a coil according to its wire cross-section, A, and length, l, is obtained from (17): where, ρ c is the copper resistivity. Moreover, the length of the wire increases due to the bunching process of Litz wire, which has a direct impact on the amount of resistance.
To consider twist factor, coefficient of β is assumed, which is considered about 1.02. Hence, the final length for the Litz wire is calculated by (18): The real cross-sectional area of the Litz wire is obtained from (19): where n is the number of strands and d st is the diameter of each strand. It is worth noting that the proposed algorithm can be used for different coils, e.g. square-types or D-types, by using their corresponding equations similar to the wheeler's equation for spiral coils, i.e. equation (15).

C. MODELING OF AC RESISTANCE
Since the switching frequency in IPT system is usually high, the skin effect and proximity effect should be considered. Therefore, it is necessary to select diameter of the strands smaller than twice the skin depth, δ, i.e. d st ≤ 2δ [28]. The number of strands is chosen based on the amount of the current. The alternating magnetic field induces eddy currents in closely adjacent conductors (as the turns of a coil). Hence, it can significantly increase the AC resistance of adjacent turns when compared to its DC resistance and this phenomenon is known as the proximity effect. To study the eddy-effects in the spiral coils, 2D FEM simulations with cylindrical symmetry are done using Maxwell simulator. Fig. 3 shows current density of the first and second turns of a flat spiral coil with 37 strands. As shown in Fig. 3 (b), the current density is non-uniform. Hence, AC resistance, R ac , value is greater than the R dc , which is modeled with R ac to R dc ratio, α, as follows: To analyze the skin and proximity effects, Poisson equation of magnetic vector potential is considered based on the following governing equation: where, A is Root Mean Squared (RMS) value of the magnetic vector potential, µ is the magnetic permeability, J s and J ind are RMS values of the source current density and induced current density, respectively. J ind , is obtained from (22): Equation (21) is rewritten as follows, using (22): where, σ is the electrical conductivity. Phasor form of equation (23) is obtained as follows: The magnetic vector potential is only in the ϕ-direction and the derivative of the magnetic vector potential in the ϕ-direction, A ϕ , with respect to ϕ-direction is zero. Therefore, the differential equation form of (24) is achieved in cylindrical coordinates as follows: After solving the magnetic vector potential, the losses should be post-processed. The losses are calculated in the natural angular frequency, ω o , as the AC losses. In addition, losses are calculated for DC currents as the DC losses. Hence, by dividing the AC losses by the DC losses, α is achieved as follows: In this paper, in order to obtain α, the eddy current analysis have been done in the Maxwell software for f o = 85 kHz. The proximity and skin effects depend on several variables, such as space between adjacent turns, S, strand diameter, d st , inner diameter of spiral coil, D i , the number of strands in Litz wire, n, and turn number of spiral coil, N . In this paper, it is approximated that these parameters are independent of each other and α is considered as function of separate variables, as derived by (27). The following assessments show that the proposed assumption has more than 95% accuracy.
where, G is correction factor, α s is AC to DC resistance ratio by changing S, α D i is the ratio by changing D i , α n is the ratio by changing n and α N is the ratio by changing N . Fig. 4 (a) shows the Maxwell results of α based on the space variations between adjacent turns from 0.5 mm up to 3 mm, while d st = 0.2 mm, D i = 200 mm, N = 20 and n = 37. According to the simulation result, a third-order interpolation has been taken as expressed in equation (28).
This procedure can be considered for each SWG strand. In this paper, the model is derived for diameters of 0.2, 0.3 and 0.35 mm, as expressed in Tables 1 up to 3, for specific  ranges included in Table 4. The ranges presented in Table 4 satisfy a wide range of practical applications. To correct the obtained equations, the correction factor, G, is achieved as 0.6678, 0.1287 and 0.0372 for strand diameters of 0.2, 0.3 and 0.35 mm, respectively.   Table 5. In general, the accuracy of the estimation equation is almost 95 %, which makes it a good candidate for universal   modeling of spiral coils' AC resistances for wireless power transfer systems.

III. DESIGN OPTIMIZATION OF THE IPT SYSTEM A. OPTIMIZATION ALGORITHM
In this paper, the proposed optimization algorithm is implemented in MATLAB based on Genetic Algorithm (GA), while other optimization algorithms can be utilized for the proposed flowchart. For GA, six variables are considered for degrees of freedom while for each variable 200 individuals are considered for population size and the maximum runs are limited to 100.
The main variables include, primary and secondary inductances, L p and L s , space between adjacent turns in primary and secondary coils, S p and S s , and inner diameter of primary and secondary coils, D p i and D s i . To show the concepts behind the proposed optimization algorithm, three optimizations are executed for different strands and powers. Regarding Fig. 5, the proposed optimization algorithm is presented in twelve steps as follows: Step 1: general specifications of the IPT system including V dc , V Bat , reference output power, P Ref , coupling factor and f o are defined by the user. Moreover, the equivalent load resistance is calculated from equation (5).
Step 2: The ranges of L p , L s , S p , S s , D p i and D s i are specified as presented in Table 6. Then GA generates random values for the variables which are between the lower and the upper bounds. Step 3: According to the equations (7), (16) and the generation values of step 2, M , C p and C s are calculated.
Step 4: To calculate r p and r s , it is necessary to determine the cross-section of the wires. At first, the r p and r s are neglected and P out and θ are calculated as functions of switching frequency for f o ≤ f s ≤ 2f o with steps of 10 Hz. Then, according to (8), the primary to secondary current ratio, k s , is obtained from (32) at the first point of f s , where P out is close to P Ref for θ ≥ 10 0 .
Therefore, according to the secondary current, derived by (4) and a proper RMS current density, J , of about 3 A/mm 2 , the cross-section area of secondary coil, A s , is obtaine.
Using A and k s , the cross-section area of the primary coil, A p , is calculated from (33): Then, regarding (19) and (33), the number of strands is obtained and rounded for the Tx and Rx, separately.
Step 5: Using (15), L p and L s , the turn number of primary and secondary coils, N p and N s , is obtained.
Step 6: The dc resistance of primary and secondary coils are obtained using equations (17) up to (19). Using the proposed equations (20) and (27) up to (31), the AC resistance is calculated.
Step 7: By knowing r p and r s from previous step, the vectors of output power, P out ( f s ), efficiency, η (f s ) and the angle of impedance of the IPT system, θ (f s ), are calculated for f o up to 2f o using the equations (8) up to (13). These vectors are stored in a lock-up table. The frequency band is from f o up to 2f o with 10 Hz steps.
Step 8: To guarantee ZVS, the first point of f s is determined at which the corresponding phase angle is greater than 10 0 , i.e. f * . Then, all components of P out (f s ) and η (f s ) that are before this point, i.e. f * , are deleted and the new lock-up table is generated from f * up to 2f o .
Step 9: The maximum component of the efficiency vector, η opt , is found. Also, the corresponding output power, P opt , and its frequency, f opt , are determined.
Step 10: Using the values of η opt and P opt from previous step, the value of the proposed objective function, Obj, is calculated from (34) as follow: In (34), λ = 1, and γ is considered 10 −4 . It is worth nothing that using multi-objective GA for the two objectives, gives more accurate results. However, for simplicity in presenting the proposed algorithm, the objective function is considered based on penalty coefficients, i.e., λ and γ .
Step 11: Then, the stop criteria of the optimization algorithm is checked regarding the predefined absolute tolerance or maximum number of iterations of the GA.
Step 12: The optimization results are printed based on the variables of IPT system. Table 6 and considering f o = 85 kHz, V Bat = 46 V, and V dc = 90 V, three different designs are presented in Table 7   Furthermore, the optimization results from MATLAB are verified in PSIM simulator for case 2. In this simulation, I p is lagging about 11 0 at f s ≈ 87.4 kHz and P in and P out are about 223.74 and 213.85 W, i.e. η ≈ 95.58 %, which is in good agreement with the MATLAB outputs. In addition, for case 3, R e , is about 1.71 .

Regarding boundaries presented in
Regarding [21] and the parameters of SS IPT system, bifurcation occurs when R e is derived by the following condition: Regarding Table 7, equation (35) and R e = 1.71 , the IPT system has a weak bifurcation which shows importance of the proposed method.
Moreover, for high power applications, e.g. fast chargers, and with large dimensions, bifurcation will occur. In bifurcation condition, the phase plot of the total impedance seen from the primary side, Z in , has three zero crossing points while the proper operating frequencies occurs larger than the third zero crossing. Hence, the proposed algorithm is devised such that search the optimum values for operating frequencies larger than the third zero crossing, as described in Fig. 5 and Step 8.

C. FURTHER DISCUSSION AND CONSIDERATIONS
A practical design requires iterative process to satisfy the all constraints. Here, the air-gap between the two coils is one of the most important parameters associated with a real IPT design. After optimization, using the proposed algorithm, a simple magneto-static FEM simulation must be done to find the air-gap associated with the coil dimensions. It is worth noting that the FEM simulation has significantly less meshes and simulation time, while each turn can be considered with no strands. For a given mutual, the inner diameters plays key role in determining the air-gap. Hence, as shown in Fig. 6, if the derived air-gap, h, from FEM simulation is less than the required air-gap, h r , the lower bound of inner diameter must be increase by ratio h r /h ratio and new optimization process must be run again, using Fig. 5. Otherwise, for h > h r , upper bound of outer diameter of coils must be decrease by ratio h/h r . The new boundaries for coils' inner diameters are achieved using (36) and (37).
where, NLB and NUB are the new lower and upper bounds for the inner diameters, respectively. OLB and OUB are the old lower and upper bounds for the inner diameters in previous step, respectively. This process is done iteratively to achieve the desired air-gap. Regarding Fig. 6, ε is considered as the allowable error for airgap design. Finally, depending on application, some parameters may be considered as predefined variables. However, as can be seen, a universal algorithm is proposed to be useful for all practical IPT systems. As an example, consider that the turn numbers are predefined variables. Hence, the Step 5 of the algorithm will be removed and by knowing N s and N p , L p and L s variables can be derived by (15) based on the other variables, i.e. D i , D w and S. As a result, the L p and L s variables will be removed.  Table 7.

IV. EXPERIMENTAL RESULTS
A laboratory prototype has been developed based on optimum parameters in section III based on case 2, as shown in Fig. 7. The half-bridge inverter is constructed by two IRFP250 power MOSFETs, S1 and S2 and IR2104 gate driver. Regarding P out and V Bat , R e , is about 7.85 .
The primary and secondary coils are constructed according to case 2 in Table 7. The Litz wires are constructed by 19 insulated strands with d st = 0.35 mm. For 22 turns, L p is about 88.5 µH, and 18 turns for the Rx coil, L s is 54.3 µH. Regarding M ≈ 13.77 µH, the air gap between the primary and secondary coils is achieved about 7 cm. Fig. 8 (a) shows the primary voltage and current at steady-state condition. The primary current lags the primary voltage about 10 0 , while input power is about 216 watts. Fig. 8 (b) shows the load current and voltage in steady-state condition, that the magnitude of load current and voltage is about 7.4 A and 54 V, respectively. RMS value of secondary current using integral formula is about 5.12 A. According to the equivalent load resistance, P out is about 205.6 W. Hence, the maximum efficiency of IPT system is about 95.2 %. Therefore, the simulation results of section III are verified by the experimental results.

V. CONCLUSION
This paper has presented a new universal design optimization algorithm for IPT systems based on series-series compensation topology and flat spiral coil. To consider AC equivalent resistance of the coils, a new approximation method using separation of variables is proposed and developed based on eddy current simulations. In the proposed algorithm, the parameters of the spiral coils are determined to achieve maximum efficiency for a desired output power, dc-link voltages and natural frequency. Moreover, the algorithm searches the maximum efficiency in the operating frequencies where ZVS occurs even under bifurcation conditions. This paper shows the importance of considering the bifurcation conditions in the optimum design algorithm, especially for high power designs. He is currently an Associate Professor of electrical engineering with Qatar University. His research interests include control systems theory, neural networks, fuzzy control, and electric drive systems. VOLUME 8, 2020