Fast Prediction of Characteristics in Wound Rotor Synchronous Condenser Using Subdomain Modeling

Abstract

Wound rotor synchronous condensers (WRSCs) are DC-excited rotor machines that utilize rotor winding instead of permanent magnets. Their voltage regulator controls the rotor field to generate or absorb reactive power, thereby regulating grid voltage or improving power factor. A key characteristic of a WRSC is the compounding curve, which shows the required rotor current under specific stator current and voltage conditions. This paper presents an approach for quickly calculating the electromagnetic parameters of a WRSC using a mathematical method. After determining magnetic flux density, induced voltage, and inductance through analytical methods, the Park and Clarke transformations are applied to derive the dq-frame quantities, enabling prediction of active and reactive powers and compounding curve characteristics. The 60 Hz model was evaluated through comparison with finite element method (FEM) simulations. Results of flux density, induced voltage, and the compounding curve under varying rotor and stator current conditions showed that the proposed method achieved comparable performance to FEM simulation while reducing computational time by half.

1. Introduction

Growing environmental concerns and the impacts of climate change are driving the shift from fossil-fuel-based power generation to renewable energy sources. Integrating renewable energy into power systems worldwide helps reduce greenhouse gas emissions and addresses the significant CO₂ emissions of the energy sector, offering sustainable and eco-friendly electricity to power electronic systems. Achieving net-zero emissions by 2050 is essential to limiting global temperature rise, prompting significant changes in the global energy landscape. By 2050, variable renewable energy sources like wind and solar power are expected to dominate the energy mix, with global electricity production projected to increase by 70%. However, the intermittency of renewable energy requires innovative solutions to maintain grid stability, including advanced energy storage and flexible demand management. Consequently, increasing the installed capacity of renewable energy sources in power systems has become a global priority [1,2,3,4].

Most renewable energy resources (RERs) connect to power systems through power electronic converters, which lack the inertia provided by traditional fossil fuel generators. The high penetration of RERs, coupled with the retirement of aging fossil-fuel-based generators, has reduced power system inertia. Moreover, most RERs are located in rural areas and connect to power systems via long transmission lines. These long lines have high impedance, which lowers the power system’s short circuit ratio [5]. As a result, transmitting active power through the system consumes significant reactive power, reducing the AC grid’s dynamic reactive power reserve and making voltage stability a critical issue [6].

A synchronous condenser is a machine used as a reactive power compensation device to enhance AC power system performance. It provides two primary functions: load compensation and voltage support. Load compensation improves power factor, balances loads, and eliminates current harmonics in nonlinear industrial loads. Voltage support aims to reduce voltage fluctuations at transmission line terminals [7,8], increase the maximum transmittable active power, and enhance AC system stability [9].

Studies in [10,11] discuss converting existing or retired synchronous generators into synchronous condensers. This conversion involves evaluating and modifying auxiliary equipment, accelerating the generator slightly above synchronous speed, powering the excitation system, and synchronizing the condenser with the transmission system. The system maintains the desired terminal voltage of the condenser using automatic voltage regulation controls.

In analyzing electrical machines, particularly synchronous condensers, air-gap flux density is crucial for designers, as it enables determination of parameters such as voltage, torque, and power. Previous studies have mainly relied on finite element method (FEM) simulations, which can be time-consuming, or on magnetic equivalent circuit approaches [12], which often lack precision. Analytical methods offer faster computation with results that align closely with FEM simulations, making them an attractive option for pre-design, analysis, and optimization of wound rotor synchronous condensers (WRSCs). Two main analytical approaches have been explored: subdomain modeling and harmonic modeling.

Subdomain modeling involves solving Maxwell’s equations and boundary/interface conditions using Fourier series expansion. This approach has been applied to surface permanent magnet (SPM) motors [13]. The second approach, harmonic modeling, addresses permeability issues using a complex Fourier series and a Cauchy product, allowing consideration of nonlinear core characteristics. It was initially introduced in analyzing switched reluctance motors (SRMs) [14]. Based on these analytical approaches, ref. [15] optimized cogging torque by investigating the slot opening and magnet pitch ratio, while [16,17,18] optimized output power in magnetically geared machines (MG) and SPMs. For multiphysics analysis, ref. [19] applied the subdomain method to predict stress and deformation in MGMs, and [20] solved for thermal distribution in the stator slot of SPMs.

To the best of our knowledge, analytical methods have not been applied to WRSCs. Thus, this study extensively analyzes a WRSC using a subdomain approach to predict air-gap flux density, voltage, and power-factor-related characteristics. Section 2 presents the main content in four subsections. The Section 2.1 introduces the application of partial differential equations (PDEs) to WRSCs. The Section 2.2 establishes an equation system to solve for unknown coefficients based on boundary conditions. The Section 2.3 defines WRSC parameters, and the fourth presents the compounding curve calculation. Finally, Section 3 presents FEM simulation results validating the proposed method.

2. Subdomain Modeling

2.1. Governing Partial Differential Equations (PDEs)

Figure 1 describes the WRSC and its simplified models employed in the subdomain modeling. Initially, the following assumptions were made:

• The end effects are ignored;
• The problem is two-dimensional (2-D) in polar coordinates;
• Magnetic vector potential A, current density J, and magnetic flux density vector B have the following non-zero components: $A=\left[\mathrm{0,0},A_z\right]$; $J=\left[\mathrm{0,0},J_z\right]$; $B=\left[B_r,B_{\theta },0\right]$;
• The core materials have infinite permeability.

The WRSC model, as illustrated in Figure 1b, is divided into five subdomains: rotor, rotor opening, air gap, stator opening, and stator slots. Each subdomain is represented by a vector potential of the form $A_z^p$, $A_z^q$, $A_z^I$, $A_z^i$, $A_z^j$, with respective harmonic components $h$, $g$, $n$, $k$, $m$. The model includes $N_r$ rotor and $N_s$ stator slots, where indices $p$, $q$, $i$, and $j$ denote specific rotor and stator slots with initial positions ${\theta }_p$, ${\theta }_q$, ${\theta }_i$, and ${\theta }_j$. A 2-D analysis in polar coordinates is used to derive partial differential equations (PDEs) for each region:

$${\displaystyle \frac{{\partial }^2{A}_z^p}{\partial r^2}}+{\displaystyle \frac{\partial A_z^p}{r\partial r}}+{\displaystyle \frac{{\partial }^2{A}_z^p}{r^2\partial {\theta }^2}}=-{\mu }_0{J}_z^p$$

(1)

$${\displaystyle \frac{{\partial }^2{A}_z^q}{\partial r^2}}+{\displaystyle \frac{\partial A_z^q}{r\partial r}}+{\displaystyle \frac{{\partial }^2{A}_z^q}{r^2\partial {\theta }^2}}=0$$

(2)

$${\displaystyle \frac{{\partial }^2{A}_z^I}{\partial r^2}}+{\displaystyle \frac{\partial A_z^I}{r\partial r}}+{\displaystyle \frac{{\partial }^2{A}_z^I}{r^2\partial {\theta }^2}}=0$$

(3)

$${\displaystyle \frac{{\partial }^2{A}_z^i}{\partial r^2}}+{\displaystyle \frac{\partial A_z^i}{r\partial r}}+{\displaystyle \frac{{\partial }^2{A}_z^i}{r^2\partial {\theta }^2}}=0$$

(4)

$${\displaystyle \frac{{\partial }^2{A}_z^j}{\partial r^2}}+{\displaystyle \frac{\partial A_z^j}{r\partial r}}+{\displaystyle \frac{{\partial }^2{A}_z^j}{r^2\partial {\theta }^2}}=-{\mu }_0{J}_z^j$$

(5)

The air gap and slot-opening regions can be modeled using Laplace’s equation, while the rotor-stator slot subdomains can be represented by Poisson’s equation. To facilitate the Fourier series expansion of the vector potential, the right-hand side of Poisson’s equation must be reformulated in Fourier form.

In the double-layer stator winding, the current density in the j-th slot is expressed as

$$\begin{gathered} {\mathit{J}}_{\mathit{z}}^{\mathit{j}}=\left(J_0^j+\sum _{m=1,2}^{\infty }{J}_m^j\mathrm{c}\mathrm{o}\mathrm{s}\left({\displaystyle \frac{m\pi }{\delta }}\left(\theta -{\theta }_j\right)\right)\right){\mathit{i}}_{\mathit{z}} \\ \begin{array}{l}\mathrm{where}\\ J_0^j={\displaystyle \frac{J_0^{j-layer1}+J_0^{j-layer2}}{2}}\\ J_m^j=\left[\begin{array}{l}{\displaystyle \frac{2}{m\pi }}\left(J_0^{j\_layer1}-J_0^{j\_layer2}\right)\sin \left({\displaystyle \frac{m\pi }{2}}\right)\leftrightarrow \mathrm{n}\mathrm{o}\mathrm{n} \mathrm{o}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{l}\mathrm{a}\mathrm{p}\mathrm{p}\mathrm{i}\mathrm{n}\mathrm{g}\\ 0 \leftrightarrow \mathrm{o}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{l}\mathrm{a}\mathrm{p}\mathrm{p}\mathrm{i}\mathrm{n}\mathrm{g}\end{array}\right.\end{array} \end{gathered}$$

(6)

In the single-layer rotor winding, the current sheet is simply described as ${\mathit{J}}_{\mathit{z}}^{\mathit{p}}=J_0^p{\mathit{i}}_{\mathit{z}}$.

2.2. General Solutions

The infinite permeability assumption in the core leads to tangential-direction boundary conditions.

$${\left.{\displaystyle \frac{\partial A_z^p}{\partial \theta }}\right|}_{\theta ={\theta }_p}={\left.{\displaystyle \frac{\partial A_z^p}{\partial \theta }}\right|}_{\theta ={\theta }_p+\alpha }$$

(7)

$${\left.{\displaystyle \frac{\partial A_z^q}{\partial \theta }}\right|}_{\theta ={\theta }_q}={\left.{\displaystyle \frac{\partial A_z^q}{\partial \theta }}\right|}_{\theta ={\theta }_q+\gamma }$$

(8)

$${\left.{\displaystyle \frac{\partial A_z^i}{\partial \theta }}\right|}_{\theta ={\theta }_i}={\left.{\displaystyle \frac{\partial A_z^i}{\partial \theta }}\right|}_{\theta ={\theta }_i+\beta }$$

(9)

$${\left.{\displaystyle \frac{\partial A_z^j}{\partial \theta }}\right|}_{\theta ={\theta }_j}={\left.{\displaystyle \frac{\partial A_z^j}{\partial \theta }}\right|}_{\theta ={\theta }_j+\delta }$$

(10)

Applying these four boundary conditions to Equations (1)–(5), general solutions of PDEs are expressed as a sum of homogeneous and particular solutions.

$$A_z^p={\mathit{A}}_{\mathbf{0}}^{\mathit{p}}+\ln \left(\mathrm{r}\right) {\mathit{B}}_{\mathbf{0}}^{\mathit{p}}-{\displaystyle \frac{{\mu }_0}{4}}{r}^2{J}_0^p+\sum _{h=1,2}^{\infty }\left(r^{h{\displaystyle \frac{\pi }{\alpha }}}{\mathit{C}}_{\mathit{h}}^{\mathit{p}}+r^{-h{\displaystyle \frac{\pi }{\alpha }}}{\mathit{D}}_{\mathit{h}}^{\mathit{p}}\right)\cos \left(h{\displaystyle \frac{\pi }{\alpha }}\left(\theta -{\theta }_p\right)\right)$$

(11)

$$A_z^q={\mathit{A}}_{\mathbf{0}}^{\mathit{q}}+\ln \left(\mathrm{r}\right) {\mathit{B}}_{\mathbf{0}}^{\mathit{q}}+\sum _{g=1,2}^{\infty }\left(r^{g{\displaystyle \frac{\pi }{\gamma }}}{\mathit{C}}_{\mathit{g}}^{\mathit{q}}+r^{-g{\displaystyle \frac{\pi }{\gamma }}}{\mathit{D}}_{\mathit{g}}^{\mathit{q}}\right)\cos \left(g{\displaystyle \frac{\pi }{\gamma }}\left(\theta -{\theta }_q\right)\right)$$

(12)

$$A_z^I=\sum _{n=1,2}^{\infty }\left(r^n{\mathit{A}}_{\mathit{n}}^{\mathit{I}}+r^{-n}{\mathit{B}}_{\mathit{n}}^{\mathit{I}}\right)\sin \left(n\theta \right)+\left(r^n{\mathit{C}}_{\mathit{n}}^{\mathit{I}}+r^{-n}{\mathit{D}}_{\mathit{n}}^{\mathit{I}}\right)\cos \left(n\theta \right)$$

(13)

$$A_z^i={\mathit{A}}_{\mathbf{0}}^{\mathit{i}}+\ln \left(\mathrm{r}\right) {\mathit{B}}_{\mathbf{0}}^{\mathit{i}}+\sum _{k=1,2}^{\infty }\left(r^{k{\displaystyle \frac{\pi }{\beta }}}{\mathit{C}}_{\mathit{k}}^{\mathit{i}}+r^{-k{\displaystyle \frac{\pi }{\beta }}}{\mathit{D}}_{\mathit{k}}^{\mathit{i}}\right)\cos \left(k{\displaystyle \frac{\pi }{\beta }}\left(\theta -{\theta }_i\right)\right)$$

(14)

$$\begin{gathered} A_z^j={\mathit{A}}_{\mathbf{0}}^{\mathit{j}}+\ln \left(\mathrm{r}\right) {\mathit{B}}_{\mathbf{0}}^{\mathit{j}}-{\displaystyle \frac{{\mu }_0}{4}}{r}^2{J}_0^j+ \\ \sum _{m=1,2}^{\infty }\left(r^{m{\displaystyle \frac{\pi }{\delta }}}{\mathit{C}}_{\mathit{m}}^{\mathit{j}}+r^{-m{\displaystyle \frac{\pi }{\delta }}}{\mathit{D}}_{\mathit{m}}^{\mathit{j}}+{\displaystyle \frac{{\mu }_0{r}^2{J}_m^j}{{\left(m\frac{\pi }{\delta }\right)}^2-4}}\right)\cos \left(m{\displaystyle \frac{\pi }{\delta }}\left(\theta -{\theta }_j\right)\right) \end{gathered}$$

(15)

In total, there are 20 coefficients in (11)–(15). Thus, to determine a unique solution requires us to derive 20 corresponding equations. The continuity of the vector potential radial components results in the following boundary conditions:

$$r=R_1\to {\displaystyle \frac{\partial A_z^p}{\partial r}}=0$$

(16)

$$r=R_2\to \left\{\begin{array}{c}{\displaystyle \frac{\partial A_z^p}{\partial r}}={\displaystyle \frac{\partial A_z^q}{\partial r}}\\ A_z^q=A_z^p\end{array}\right.$$

(17)

$$r=R_3\to \left\{\begin{array}{c}{A}_z^q=A_z^I\\ {\displaystyle \frac{\partial A_z^I}{\partial r}}={\displaystyle \sum _{p=1,2}^{N_r}}{\displaystyle \frac{\partial A_z^p}{\partial r}}\end{array}\right.$$

(18)

$$r=R_4\to \left\{\begin{array}{c}{\displaystyle \frac{\partial A_z^I}{\partial r}}={\displaystyle \sum _{i=1,2}^{N_s}}{\displaystyle \frac{\partial A_z^i}{\partial r}}\\ A_z^i=A_z^I\end{array}\right.$$

(19)

$$r=R_5\to \left\{\begin{array}{c}{A}_z^i=A_z^j\\ {\displaystyle \frac{\partial A_z^j}{\partial r}}={\displaystyle \frac{\partial A_z^i}{\partial r}}\end{array}\right.$$

(20)

$$r=R_6\to {\displaystyle \frac{\partial A_z^j}{\partial r}}=0$$

(21)

By employing integrals $sni(n,{\theta }_i,\beta )$, $rni(n,{\theta }_i,\beta )$, $gkni\left(k,n,{\theta }_i,\beta \right)$, $fkni\left(k,n,{\theta }_i,\beta \right)$, and $Fmk\left(m,k,\beta ,\delta \right)$ in [19] to shorten equations, 20 equations are given in Appendix A and Appendix B.

2.3. Matrix Representation and Solving Electromagnetic Quantities

Reformulating the above equations into matrix and vector forms makes it possible to obtain analytical solutions using MATLAB R2022b. For instance, in a WRSM model with $N_r=32$ rotor and $N_s=42$ stator slots and harmonic order numbers $G=H=K=M=5$ and $N=100$, the result is a column matrix X with $4\left(N+N_{r}G+N_{r}H+N_{s}K+N_{s}M\right)=3360$ elements.

$$\mathit{X}={\left[{\mathit{A}}_{\mathit{n}}^{\mathit{I}}{\mathit{B}}_{\mathit{n}}^{\mathit{I}}{\mathit{C}}_{\mathit{n}}^{\mathit{I}}{\mathit{D}}_{\mathit{n}}^{\mathit{I}}{\mathit{A}}_{\mathbf{0}}^{\mathit{p}}{\mathit{B}}_{\mathbf{0}}^{\mathit{p}}{\mathit{C}}_{\mathit{h}}^{\mathit{p}}{\mathit{D}}_{\mathit{h}}^{\mathit{p}}{\mathit{A}}_{\mathbf{0}}^{\mathit{q}}{\mathit{B}}_{\mathbf{0}}^{\mathit{q}}{\mathit{C}}_{\mathit{g}}^{\mathit{q}}{\mathit{D}}_{\mathit{g}}^{\mathit{q}}{\mathit{A}}_{\mathbf{0}}^{\mathit{i}}{\mathit{B}}_{\mathbf{0}}^{\mathit{i}}{\mathit{C}}_{\mathit{k}}^{\mathit{i}}{\mathit{D}}_{\mathit{k}}^{\mathit{i}}{\mathit{A}}_{\mathbf{0}}^{\mathit{j}}{\mathit{B}}_{\mathbf{0}}^{\mathit{j}}{\mathit{C}}_{\mathit{m}}^{\mathit{j}}{\mathit{D}}_{\mathit{m}}^{\mathit{j}}\right]}^{\mathit{T}}$$

(22)

After solving the coefficients, flux density at the air gap (radius is $R_e$) is derived as:

$$B_{\theta }=-{\displaystyle \frac{\partial A^I}{\partial r}}=-{\displaystyle \frac{1}{R_e}}\sum _{n=1,2}^{\infty }\left({R_e}^n{\mathit{A}}_{\mathit{n}}^{\mathit{I}}-{R_e}^{-n}{\mathit{B}}_{\mathit{n}}^{\mathit{I}}\right)\mathit{sin}\left(n\theta \right)n+\left({R_e}^n{\mathit{C}}_{\mathit{n}}^{\mathit{I}}-{R_e}^{-n}{\mathit{D}}_{\mathit{n}}^{\mathit{I}}\right)\mathit{cos}\left(n\theta \right)n$$

(23)

$$B_r={\displaystyle \frac{1}{R_e}}{\displaystyle \frac{\partial A^I}{\partial \theta }}={\displaystyle \frac{1}{R_e}}\sum _{n=1,2}^{\infty }\left({R_e}^n{\mathit{A}}_{\mathit{n}}^{\mathit{I}}+{R_e}^{-n}{\mathit{B}}_{\mathit{n}}^{\mathit{I}}\right)\cos \left(n\theta \right)n-\left({R_e}^n{\mathit{C}}_{\mathit{n}}^{\mathit{I}}+{R_e}^{-n}{\mathit{D}}_{\mathit{n}}^{\mathit{I}}\right)\sin \left(n\theta \right)n$$

(24)

To compute the induced voltage in a three-phase motor with stack length $L_{stk}$, the flux through each slot cross-section at a given rotor position ${\theta }_0$ was calculated. Uniform current density across the slot area was assumed, allowing for the flux in the j-th slot to be determined by integrating the vector potential across the slot area, as given in Equation (15).

$$\begin{gathered} {\mathsf{\Phi }}_j=T_{pole}{\displaystyle \frac{l_{stk}}{A_{slot}}}\iint A_r^j\left(\theta ,z\right)rdrd\theta \\ ={\displaystyle \frac{L_{stk}}{\frac{R_6^2-R_5^2}{2}}}\left({\mathit{A}}_{\mathbf{0}}^{\mathit{j}}{\displaystyle \frac{R_6^2-R_5^2}{2}}+{\mathit{B}}_{\mathbf{0}}^{\mathit{j}}\left(\ln \left(R_6\right){\displaystyle \frac{R_6^2}{2}}-\ln \left(R_5\right){\displaystyle \frac{R_5^2}{2}}+{\displaystyle \frac{R_6^2-R_5^2}{4}}\right)-J_0^j{\displaystyle \frac{{\mu }_0\left(R_6^4-R_5^4\right)}{16}}\right) \end{gathered}$$

(25)

The phase flux is calculated as the sum of fluxes in slots associated with each phase. A connecting matrix $\left[C\right]$ with the dimensions ${3\times N}_s$ was used to present the stator winding distribution in the slots, as shown in [13], where indices 1 and −1 denote positive and negative, and 0 denotes a phase absence in the slot. $T_{slot}$ represents the phase-winding turns per slot, and the phase fluxes are expressed as:

$${\left[\begin{array}{c}{\Phi }_A\\ {\Phi }_B\\ {\Phi }_C\end{array}\right]}_{3\times 1}=T_{slot}{\left[C\right]}_{{3\times N}_s}{\left[{\Phi }_j\right]}_{N_s\times 1}$$

(26)

The induced voltage can be defined at a given rotor speed ${\omega }_m$ as follows:

$${\left[\begin{array}{c}{U}_A\\ U_B\\ U_C\end{array}\right]}_{3\times 1}={\displaystyle \frac{N_r}{2}}{\omega }_m{\displaystyle \frac{d}{d{\theta }_0}}{\left[\begin{array}{c}{\Phi }_A\\ {\Phi }_B\\ {\Phi }_C\end{array}\right]}_{3\times 1}$$

(27)

Notably, without stator current contribution in the no-load condition, the induced voltage becomes BEMF.

2.4. Compounding Curve

Figure 2 illustrates three compounding curves, where the red, blue, and orange lines represent power factors of 1 and 0.95 leading (WRSC generates reactive power), and 0.95 lagging (WRSC absorbs reactive power), respectively. For instance, at a stator load of 105 A, the rotor generates 86 A to achieve a power factor of 1. If the rotor current exceeds 86 A, the machine becomes overexcited, and excess reactive power flows to the grid. Conversely, if the rotor current is below 86 A, the machine becomes under-excited, lowering its operating voltage below the grid voltage, prompting the grid to supply reactive power to compensate for the power deficit in the machine.

The compounding curve is essential for assessing WRSC performance, allowing de-signers to estimate the field current necessary for a given stator load to achieve a desired electrical system power factor. This section presents the procedure for deriving the unity power factor line and applies it to salient and non-salient rotors. The procedure involves:

• Transforming phase quantities (flux, voltage, and current) obtained from subdomain modeling into the dq-frame under current excitation.
• Deriving and transforming self- and mutual inductances into the dq-frame.
• Building and solving a d-axis current function with a reactive power of zero from the voltage and power equations in the dq-frame.

By aligning the rotor pole with phase A in the initial position, the quantities of voltage, current, flux, and inductance can be transferred to the dq-frame using Park and Clarke transformations [21]. By assuming that machine saturation is neglected, Park and Clarke transformations are given as follows:

$${UI\lambda }_{dqz}={\displaystyle \frac{2}{3}}\left[\begin{array}{ccc}\mathrm{c}\mathrm{o}\mathrm{s}\left({\theta }_0\right)& \cos \left({\theta }_0-{\displaystyle \frac{2\pi }{3}}\right)& \cos \left({\theta }_0+{\displaystyle \frac{2\pi }{3}}\right)\\ -\mathrm{s}\mathrm{i}\mathrm{n}\left({\theta }_0\right)& -\sin \left({\theta }_0-{\displaystyle \frac{2\pi }{3}}\right)& -\sin \left({\theta }_0+{\displaystyle \frac{2\pi }{3}}\right)\\ {\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}& {\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}& {\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\end{array}\right]\left[\begin{array}{c}{UI\lambda }_a\\ {UI\lambda }_b\\ {UI\lambda }_c\end{array}\right]$$

(28)

$$\left[\begin{array}{ccc}{L}_d& M_{dq}& M_{dz}\\ M_{qd}& L_q& M_{qz}\\ M_{zd}& M_{zq}& L_z\end{array}\right]=\sqrt{{\displaystyle \frac{2}{3}}}\left[\begin{array}{ccc}1& \cos \left({\displaystyle \frac{2\pi }{3}}\right)& \cos \left({\displaystyle \frac{4\pi }{3}}\right)\\ 0& \sin \left({\displaystyle \frac{2\pi }{3}}\right)& \sin \left({\displaystyle \frac{4\pi }{3}}\right)\\ {\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\sqrt{2}$}\right.}& {\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\sqrt{2}$}\right.}& {\displaystyle \raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\sqrt{2}$}\right.}\end{array}\right]\left[\begin{array}{ccc}{L}_{aa}& M_{ab}& M_{ac}\\ M_{ba}& L_{bb}& M_{bc}\\ M_{ca}& M_{cb}& L_{cc}\end{array}\right]\sqrt{{\displaystyle \frac{2}{3}}}$$

(29)

To calculate the winding self-inductance and mutual inductances in (29), the authors apply an excitation current of 1 A to the first winding, while setting the currents in all other windings to zero, excluding rotor current effects from the analysis. For example, phase A is set to 1 A, while phases B and C and the rotor current are set to zero. The obtained flux linkages of phases A, B, and C (${\psi }_a,{\psi }_b,{\psi }_c$) are then used to calculate the inductance as follows:

$$\left[L_{aa} L_{ab} L_{ac}\right]={\displaystyle \frac{\left[{\psi }_a {\psi }_b {\psi }_c\right]}{1A}}$$

(30)

Additionally, rotor flux linkage (${\psi }_f$), used to derive rotor inductance, $M_{df}={\psi }_f/I_f$, is calculated by setting all stator winding currents to zero and applying only a 1 A current to the rotor. The flux linkage of phase A is considered the rotor flux linkage.

In high-power machines, reactance is dominant, and stator resistance can be neglected; therefore, the voltage equation is given by:

$$\left[\begin{array}{c}{U}_d\\ U_q\end{array}\right]\approx {\omega }_e\left[\begin{array}{c}-L_q{i}_q\\ L_d{i}_d+M_{df}{i}_f\end{array}\right]$$

(31)

Consequently, the active and reactive powers are defined as:

$$\begin{gathered} P=1.5(U_d{i}_d+U_q{i}_q) \\ Q=1.5(U_q{i}_d-U_d{i}_q) \end{gathered}$$

(32)

Table 1. Machine mode definition.

Mode	Power Factor	P	Q
Motor	Lagging	>0	>0
Motor	Leading	>0	<0
Generator	Lagging	<0	>0
Generator	Leading	<0	<0

shows the operating modes of the machine’s power. To operate as a capacitor bank, the machine generates reactive power (negative Q). This is achieved by applying a sufficiently large negative d-axis current along with a specific rotor current. For a constant stator current, this capacitor-like operation can also be achieved by applying a sufficient rotor current.

To achieve a unity power factor (PF = 1), commutation in WRSC is defined by a zero reactive power, Q = 0. This condition, combined with the voltage constraint $v_d^2+v_q^2=V_{smax}^2$, allows the d-axis stator current to be modeled as a function of rotor current.

$$\left(L_d^2-L_d{L}_q\right)i_d^2+M_{df}{i}_f\left(2L_d-L_q\right)i_d+{M_{df}{i}_f}^2-{\displaystyle \frac{V_{smax}^2}{{\omega }_e^2}}=0$$

(33)

By adjusting the rotor current $i_f$, the d-axis current can be solved using Equation (33), and the q-axis current can be obtained from Equation (32). These computations form a compounding curve that maps the current relationship between the rotor and stator. In salient-pole rotor machines, where dq-frame inductances are nearly identical, Equation (33) reduces to a first-order equation, simplifying analysis.

3. FEM Simulation Comparison

To validate these principles, specifications for a one-pole-pair WRSC, as listed in

Table 2. WRSC parameters.

Quantity	Symbol	Unit	Value
Inner rotor slot radius	$R_1$	mm	269.50
Outer rotor radius	$R_2$	mm	425.25
Inner stator slot radius	$R_3$	mm	525.25
Outer stator slot radius	$R_4$	mm	700.50
Outer stator radius	$R_5$	mm	1119.50
Stack length	$L_{stk}$	mm	6402.40
Vacuum permeability	${\mu }_0$	$\mathrm{k}\mathrm{g}·\mathrm{m}·{\mathrm{s}}^{-2}·{\mathrm{A}}^{-2}$	$4\pi \times {10}^{-7}$
Rotor slot number	$N_r$	-	32
Stator slot number	$N_s$	-	42
Rotor slot pitch ratio	$\alpha$	$\pi /180$ rad	4.64
Rotor opening slot pitch ratio	$\gamma$	$\pi /180$ rad	4.64
Stator opening slot pitch ratio	$\beta$	$2\pi /N_s$ rad	0.46
Stator slot pitch ratio	$\delta$	$2\pi /N_s$ rad	0.46
Rotor slot winding turn	$T_{slot}^r$	-	5/7 *
Stator slot winding turn	$T_{slot}^s$	-	2

* The 1st, 16th, 17th, and 32nd slots each contain five turns. The remaining slots accommodate seven turns.

, were simulated using a 2-D finite element method (FEM) model with the mesh setting depicted in Figure 3a. Due to mechanical constraints in MW-class WRSC, the air gap was intentionally designed to be large (up to 100 mm). The FEM simulation took 66 s to compute flux density, flux at the air gap, and induced voltage (see Figure 3b and Figure 4). By contrast, the subdomain modeling approach achieved these calculations in 33 s, showcasing the computational efficiency of the proposed method in WRSC pre-design.

Figure 3b illustrates the flux density distribution at current conditions $i_f=1 \mathrm{k}\mathrm{A};i_d=i_q=-0.5 \mathrm{k}\mathrm{A}$. The salient-pole rotor caused flux linkage to concentrate in the rotor pole and teeth near the pole, inducing mild saturation in the teeth. Under typical conditions, machines avoid saturation, making the assumption of infinite permeability in subdomain modeling suitable. However, when the field current or stator current is increased, saturation intensifies, resulting in accuracy discrepancies in the subdomain method (see Figure 4e,f).

The compounding curves derived from the subdomain method were closely aligned with the FEM results, as seen in Figure 5. To maintain a unity power factor while increasing stator current, the field current must also rise. Additionally, rotor current must increase at higher stator voltages. Calculations for these curves, derived from Equation (33), involve computing the dq-axis currents for each rotor current, then determining the stator current. The analysis covered phase peak voltages of 10, 15, 18, and 22 kV, revealing that the machine can support stator currents up to 5 kA with a maximum field current of 1.3 kA. The mathematical method obtained these results in under 1 s, which represents a significant time advantage in deriving compounding curves.

4. Conclusions

Overall, the proposed mathematical method offers a fast and effective approach for characterizing WRSC properties, including flux density, induced voltage, inductance, and compounding curves. Under low rotor and stator current conditions, the subdomain method aligns well with 2-D FEM results. However, as rotor and stator currents increase, saturation can induce errors of up to 10% in the flux density and induced voltage results. The unity-power-factor compounding curves were accurate across varying conditions. These results were obtained because of the unity power factor and suppression of demagnetization by the d-axis stator current, which weakens flux linkage and limits saturation.

Compared to 2-D FEM simulations, the mathematical method cut computation time from 66 to 33 s, facilitating faster optimizations in the WRSC design process. The reduced time can lead to exponential efficiency gains, enabling designers to employ optimization techniques more effectively.

Future research based on these findings could explore:

• Nonlinear magnetic material characteristics;
• Deriving leading and lagging compounding curves for power factor;
• Optimization techniques for a more efficient WRSC design;
• Experimental verification of the subdomain method on a test bench.