*2.7. N = 34 New Magic Number as a Consequence of the Shell Evolution*

Among various cases of shell evolution, a notable impact was made by predicting a new magic number *N* = 34. Figure 7 displays the shell evolution of some neutron orbits from Ni back to Ca isotopes, as *Z* decreases from 28 to 20. The 1 *f*5/2 orbit is between the 2*p*3/2 and 2*p*1/2 in Mayer–Jensen's shell model (see Figure 1). By loosing eight protons lying in the 1 *f*7/2 orbit of Ni isotopes (blue circles in Figure 7 left), this canonical shell structure is destroyed as the 1 *f*5/2 orbit moves up above the 2*p*1/2 orbit. This movement of 1 *f*5/2 orbit creates the *N* = 32 gap as a byproduct [44]. The energy shift in the 1 *f*5/2 orbit is due to the central and tensor forces by almost equal amounts. We mention that the *N* = 34 magic gap would not appear, if the Mayer–Jensen scheme holds, as expected, in Ni isotopes but this shift did not occur. The appearance of the *N* = 34 magic number was predicted as a result of a spin–isospin interaction in [45]. However, 12 years were required [46] until the experimental verification became feasible [47] (see Figure 7 right). The measured 2+ energy levels are included in Figure 2b. More details are presented in [15]. Further evidences have been obtained recently by different experimental probes as reported in [48,49].

**Figure 7. Left**: Schematic illustration of the shell evolution from Ni back to Ca for neutron orbits. Blue circles denote protons. The wavy line is the interaction between the proton 1 *f*7/2 orbit and the neutron 1 *f*5/2 orbit. The numbers in circles indicate magic numbers. Taken from Figure 3 of [24]. **Right**: Observed excitation energies of the 2+ 1 states. Taken from Figure 2c of [47].


Through a comprehensive analysis of AHP, the score of each solution on the Pareto frontier was acquired. The solution with the greatest score is the best one. As presented by Figures 15 and 16, the optimal design strategy is number 96 of the Pareto frontiers, with triple redundancy configuration used in the MCES. The reliability of our designed configuration increased by 5.89% and 55.34% respectively, compared with the dual redundancy and single redundancy configurations. The results indicated that our approach combining multi-objective optimization and decision-making could realize the multi-objective optimization design of an MECS. The solutions on the Pareto frontier acquired by NSGA-II arithmetic were uniformly distributed, which suggested that we could offer designers more diversified choices. Subsequently, AHP was employed for the eventual decision analyses, which enabled designers to introduce predilections and experiences for different goals.

## **5. Conclusions and Future Work**

The present research demonstrated the multi-objective optimization for an MECS. By designing the level length, the objectives are optimized. However, it is hard to minimize every goal simultaneously due to the conflicts between these objectives. Thus, an appropriate multi-objective optimization algorithm, which can compromise multiple objectives, is required to acquire the Pareto optimum solution set. Besides, since the Pareto frontier solution is not a solution but a solution set, an approach of multi-objective decision-making analysis is required to identify the optimum solution in the Pareto frontier. Herein, a combination of multi-objective optimization and multi-objective decision-making was utilized. The NSGA-II was employed to determine the optimum solution set, i.e., the Pareto frontier. Subsequently, the optimum design was acquired via an AHP. The outcome revealed the practicability of our optimization approach. In addition, the present research discovered that the multi-objective optimization and multi-objective decision-making approaches could be utilized in the optimization design of an MECS, and that this approach suited other components and systems as well.

Future work should encompass the improvement of real-time performance for an optimization algorithm, with more evaluation indicators considered.

**Author Contributions:** Conceptualization, Z.L. and R.C.; methodology, Z.L. and R.C.; software, D.L.; validation, Z.L. and D.L.; writing—original draft preparation, Z.L.; writing—review and editing, S.W.; supervision, S.W. and J.S.; project administration, J.S. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (No.<PHONE_NUMBER>0) and the program of China Scholarship Council (No.<PHONE_NUMBER>06).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** This study does not report any data.

**Conflicts of Interest:** The authors declare no conflict of interest.

## **References**

## *Article* **Novel Approach to Fault-Tolerant Control of Inter-Turn Short Circuits in Permanent Magnet Synchronous Motors for UAV Propellers**

**Aleksander Suti \* , Gianpietro Di Rito and Roberto Galatolo**

Department of Civil and Industrial Engineering, University of Pisa, Largo Lucio Lazzarino 2, 56122 Pisa, Italy<EMAIL_ADDRESS>(G.D.R<EMAIL_ADDRESS>(R.G.)

**\*** Correspondence<EMAIL_ADDRESS>Tel.: +39-05-0221-7211

**Abstract:** This paper deals with the development of a novel fault-tolerant control technique aiming at the diagnosis and accommodation of inter-turn short circuit faults in permanent magnet synchronous motors for lightweight UAV propulsion. The reference motor is driven by a four-leg converter, which can be reconfigured in case of a phase fault by enabling the control of the central point of the motor Y-connection. A crucial design point entails the development of fault detection and isolation (FDI) algorithms capable of minimizing the failure transients and avoiding the short circuit extension. The proposed fault-tolerant control is composed of two sections: the first one applies a novel FDI algorithm for short circuit faults based on the trajectory tracking of the motor current phasor in the Clarke plane; the second one implements the fault accommodation, by applying a reference frame transformation technique to the post-fault commands. The control effectiveness is assessed via nonlinear simulations by characterizing the FDI latency and the post-fault performances. The proposed technique demonstrates excellent potentialities: the FDI algorithm simultaneously detects and isolates the considered faults, even with very limited extensions, during both stationary and unsteady operating conditions. In addition, the proposed accommodation technique is very effective in minimizing the post-fault torque ripples.

**Keywords:** all-electric propulsion; electric machines; fault diagnosis; fault-tolerant control; inter-turn short circuit; modelling; simulation

## **1. Introduction**

The design of next-generation long-endurance UAVs is undoubtedly moving towards the full-electric propulsion. Though immature today, full-electric propulsion systems are expected to obtain large investments in the forthcoming years, to replace conventional internal combustion engines, as well as hybrid or hydrogen-based solutions. Full-electric propulsion would guarantee smaller CO2-emissions, lower noise, a higher efficiency, a reduced thermal signature (crucial for military applications), and simplified maintenance. However, several reliability and safety issues are still open, especially for long-endurance UAVs flying in unsegregated airspaces.

In this context, the Italian Government and the Tuscany Regional Government cofunded the project TERSA (*Tecnologie Elettriche e Radar per Sistemi aeromobili a pilotaggio remoto Autonomi*) [1], led by Sky Eye Systems (Italy) in collaboration with the University of Pisa and other Italian industries. The TERSA project aims to develop an Unmanned Aerial System (UAS) with a fixed-wing UAV, Figure 1, having the following main characteristics:

**Citation:** Suti, A.; Di Rito, G.; Galatolo, R. Novel Approach to Fault-Tolerant Control of Inter-Turn Short Circuits in Permanent Magnet Synchronous Motors for UAV Propellers. *Aerospace* **2022**, *9*, 401. https://doi.org/10.3390/ aerospace9080401

Academic Editor: Cengiz Camci

Received: 1 July 2022 Accepted: 25 July 2022 Published: 26 July 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 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/).

• Landing system: parachute and airbags; Propulsion system: Permanent Magnet Synchronous Motor (PMSM) powering a twin-blade fixed-pitch propeller;

 Take-off system: pneumatic launcher; Landing system: parachute and airbags;

• Propulsion system: Permanent Magnet Synchronous Motor (PMSM) powering a twin-blade fixed-pitch propeller; Innovative sensing systems:

*Aerospace* **2022**, *8*, x FOR PEER REVIEW 2 of 25

• Innovative sensing systems: o Synthetic aperture radar, to support surveillance missions in adverse environ-

Range: >3 km;

**Figure 1.** Rendering of the TERSA UAV. **Figure 1.** Rendering of the TERSA UAV.

With reference to the activities related to the TERSA UAV propulsion system development (which this work refers to), a special attention has been devoted to the architecture definition for obtaining fault-tolerant capabilities. In particular, propulsion systems based on three-phase PMSMs driven by conventional three-leg converters typically have failure rates from 100 to 200 per million flight hours [2], which are not compatible with the reliability and safety levels required for the airworthiness certification [3]. The failure rate of PMSMs is mostly driven by motor phase faults and converter faults (open-switch in a converter leg, open-phase, inter-turn, phase-to-leg, phase-to-ground, or capacitor shortcircuits [4]), which cover from 60% to 70% of the total fault modes [5]. Provided that the weight and envelopes required by UAV applications impede the extensive use of hardware redundancy (e.g., redundant motors), the reliability enhancement of full-electric pro-With reference to the activities related to the TERSA UAV propulsion system development (which this work refers to), a special attention has been devoted to the architecture definition for obtaining fault-tolerant capabilities. In particular, propulsion systems based on three-phase PMSMs driven by conventional three-leg converters typically have failure rates from 100 to 200 per million flight hours [2], which are not compatible with the reliability and safety levels required for the airworthiness certification [3]. The failure rate of PMSMs is mostly driven by motor phase faults and converter faults (open-switch in a converter leg, open-phase, inter-turn, phase-to-leg, phase-to-ground, or capacitor short-circuits [4]), which cover from 60% to 70% of the total fault modes [5]. Provided that the weight and envelopes required by UAV applications impede the extensive use of hardware redundancy (e.g., redundant motors), the reliability enhancement of full-electric propulsion systems can be achieved only through the redundancy of phases or stator modules [6] or by unconventional converters. By reducing design complexity, weight and envelopes, the use of four-leg converters [7–9], Figure 2, could represent a suitable solution for UAV applications.

pulsion systems can be achieved only through the redundancy of phases or stator modules [6] or by unconventional converters. By reducing design complexity, weight and envelopes, the use of four-leg converters, [7–9], Figure 2, could represent a suitable solution for UAV applications. In this converter topology, a couple of power switches are added, as stand-by devices, to the conventional three-leg bridge, enabling the control of the central point of the "star" connection (often Y-connection) of the motor. Nevertheless, to benefit from their faulttolerance capability, these devices must be promptly reconfigured when a fault occurs, and extremely fast Fault-Detection and Isolation (FDI) algorithms must be developed. Especially for PMSMs operating at high speeds (such as the ones used for aircraft propellers), phase faults can generate abnormal torque ripples and the related failure transients can potentially cause unsafe conditions [10,11].

**Figure 2.** Four-leg converter driving a three-phase PMSM with access to the central point. **Figure 2.** Four-leg converter driving a three-phase PMSM with access to the central point.

In this converter topology, a couple of power switches are added, as stand-by devices, to the conventional three-leg bridge, enabling the control of the central point of the "star" connection (often Y-connection) of the motor. Nevertheless, to benefit from their Many pieces of research on the development of FDI algorithms for PMSMs have been carried out in the last decades, and a special attention has been recently addressed to the Inter-Turn Short Circuit (ITSC) fault [12–15], which typically contributes to 10% of the PMSM failure rate [5,9,16].

fault-tolerance capability, these devices must be promptly reconfigured when a fault occurs, and extremely fast Fault-Detection and Isolation (FDI) algorithms must be developed. Especially for PMSMs operating at high speeds (such as the ones used for aircraft Faiz [12] carried out a comprehensive literature review on FDI methods for ITSC faults in PMSMs. Depending on the measurements processed by the algorithms, they can be categorized as torque-based [17,18], flux-based [19,20], electromagnetic parametersbased [21], voltage-based, and current-based.

propellers), phase faults can generate abnormal torque ripples and the related failure transients can potentially cause unsafe conditions [10,11]. Many pieces of research on the development of FDI algorithms for PMSMs have been carried out in the last decades, and a special attention has been recently addressed to the Inter-Turn Short Circuit (ITSC) fault [12–15], which typically contributes to 10% of the PMSM failure rate [5,9,16]. Faiz [12] carried out a comprehensive literature review on FDI methods for ITSC faults in PMSMs. Depending on the measurements processed by the algorithms, they can be categorized as torque-based [17,18], flux-based [19,20], electromagnetic parametersbased [21], voltage-based, and current-based. The most common strategy is given by the voltage-based methods. Immovilli and Bianchini [22] propose reconstructing harmonic patterns of the voltage components in the Clarke plane. The technique works very well in transient conditions, but it is not able to identify the shorted phase. The FDI algorithm proposed in [23], implements an estimator of back electromotive force and detects the ITSC by comparing the estimates to a reference The most common strategy is given by the voltage-based methods. Immovilli and Bianchini [22] propose reconstructing harmonic patterns of the voltage components in the Clarke plane. The technique works very well in transient conditions, but it is not able to identify the shorted phase. The FDI algorithm proposed in [23], implements an estimator of back electromotive force and detects the ITSC by comparing the estimates to a reference model. The method is very robust against load disturbances, but it requires a very accurate model of the Back Electromotive Force (BEMF). Other approaches use as fault symptoms the harmonic components of the Zero Sequence Voltage Component (ZSVC) signal. Urresty et al. [24] proposes a Vold–Kalman Filtering Order Tracking (VKF-OT) algorithm to track the ZSVCs first harmonic component. The method offers reliable indicators, which are applicable in transient phases too, but it needs to be extended in terms of fault isolation. A robust algorithm is proposed by Hang et al., in [25], where the fault indicators are defined to remove the influence of rotor speed variation and a frequency tracking algorithm based on reference frame transformation is used to extract the fault symptoms. A possible drawback is the dependency of the fault indicator on the PMSM parameters, which may vary with environmental temperature. Boileau [26] considered as a fault indicator the amplitude of the positive voltage sequence while Meinguet [27] used as an indicator the ratio between the amplitudes of negative and positive voltage sequences. These methods are not able to locate the shorted phase, neither.

model. The method is very robust against load disturbances, but it requires a very accurate Concerning the current-based methods, the most relevant one is based on the socalled Park Vector Approach (PVA) [28–31]. When the PMSM is operating in normal

model of the Back Electromotive Force (BEMF). Other approaches use as fault symptoms the harmonic components of the Zero Sequence Voltage Component (ZSVC) signal. Ur-

are applicable in transient phases too, but it needs to be extended in terms of fault isolation. A robust algorithm is proposed by Hang et al., in [25], where the fault indicators are defined to remove the influence of rotor speed variation and a frequency tracking algorithm based on reference frame transformation is used to extract the fault symptoms. A possible drawback is the dependency of the fault indicator on the PMSM parameters, which may vary with environmental temperature. Boileau [26] considered as a fault indicator the amplitude of the positive voltage sequence while Meinguet [27] used as an conditions, the current phasor in the Clarke plane draws a circular trajectory, while it follows an elliptic one if a stator fault occurs. In PVAs, the fault index is thus defined by a geometrical distortion factor, obtained by currents measurements. Though demonstrating to be excellent for online ITSC FDI, the PVA is not able to locate the shorted phase.

All the current methods have in common the observation that in case of ITSCs, the stator currents are not balanced, and higher harmonic components in the signals can be used as fault symptoms. Different signal processing methods are employed to extract these symptoms (or patterns). Frequency-domain transforms such as Fast-Fourier Transforms (FFT) are applied in [13,32–34], while mixed frequency/time-domain transforms such as Short-Time Fourier Transform (STFT) are used in [14], Wavelet Transform (WT) in [15,35,36], and Hilbert–Huang Transform (HHT) in [37,38]. Frequency-domain transforms entail the loss of transient events such as speed or loads variations, so FFT analyses, though effective in stationary conditions, must be avoided if the FDI is required during transients. By requiring a higher computational effort, STFT analyses partially compensate the FFT limits, but FDI capabilities are still limited in unsteady conditions. On the other hand, WT methods permit to decompose time domain signals into stationary and non-stationary contributions, but the appropriate choice of wavelets is a drawback. Many limits of WT techniques are removed by HHT techniques, even if they are effective when applied to signals characterised by a narrow frequency content.

Other methods, based on artificial intelligence applications (Convolutional Neural Networks are used in [39]) have been also explored, but, though the results are remarkable, these techniques are not preferred for airworthiness certification [40].

Within this research context, this paper aims to develop an innovative Fault-Tolerant Control (FTC), capable of prompt FDI and accommodation of ITSC faults in three-phase PMSMs for lightweight UAVs propulsion, in both stationary and transient conditions. The proposed FTC is composed of two sections: the first one addresses the FDI problem via an original current-based method (here referred as Advanced-PVA, APVA), the second is dedicated to the fault accommodation, obtained by activating the stand-by leg of a four-leg converter and controlling the central point of the motor Y-connection. It is worth noting that the ITSC accommodation is developed by applying the Rotor Current Frame Transformation Control (RCFTC) technique, already applied by the authors in [9] for the FTC of open-circuit faults.

The paper is articulated as follows: in the first part, the nonlinear model of the propulsion system is presented; successively, the FDI algorithms and the fault accommodation technique are described by highlighting their basic design criteria; finally, a summary of the simulation results is proposed, by characterizing the failure transients and the post-fault behaviours after the injection of ITSC faults.

## **2. Materials and Methods**

## *2.1. System Description*

The reference propulsion system, designed for the full-electric propulsion of a lightweight UAV, is composed of (Figure 3):

	- # Three-phase surface-mounted PMSM with phase windings in Y-connection;
	- # Twin-blade fixed-pitch propeller [41];
	- # Mechanical coupling joint.
	- # CONtrol/MONitoring (CON/MON) module, for the implementation of the closed-loop control and health-monitoring functions;
	- # Four-leg converter;
	- # Three current sensors (CSa, CSb, CSc), one per each motor phase;
	- # One Angular Position Sensor (APS), measuring the motor angle;

o Three current sensors (CSa, CSb, CSc), one per each motor phase; o One Angular Position Sensor (APS), measuring the motor angle;

*Aerospace* **2022**, *8*, x FOR PEER REVIEW 5 of 25

o Four-leg converter;

The CON module operates the closed-loop control of the system by implementing two nested loops on propeller speed and motor currents (via Field-Oriented Control, FOC). All the regulators implement proportional/integral actions on the tracking error signals, plus anti-windup functions with back-calculation algorithms to compensate for commands saturation. The MON module executes the health-monitoring algorithms, including the FTC proposed in this work. The CON module operates the closed-loop control of the system by implementing two nested loops on propeller speed and motor currents (via Field-Oriented Control, FOC). All the regulators implement proportional/integral actions on the tracking error signals, plus anti-windup functions with back-calculation algorithms to compensate for commands saturation. The MON module executes the health-monitoring algorithms, including the FTC proposed in this work.

### *2.2. Model of the Aero-Mechanical Section 2.2. Model of the Aero-Mechanical Section*

The dynamics of the aero-mechanical section of the propulsion system, providing the UAV with the thrust is schematically depicted in Figure 4a and is modelled by [7,9]: The dynamics of the aero-mechanical section of the propulsion system, providing the UAV with the thrust is schematically depicted in Figure 4a and is modelled by [7,9]:

$$\begin{cases} I\_p \ddot{\theta}\_p = -Q\_p - \mathbb{C}\_{gb} \left( \dot{\theta}\_p - \dot{\theta}\_m \right) - \mathbb{K}\_{gb} \left( \theta\_p - \theta\_m \right) + Q\_d\\ I\_m \ddot{\theta}\_m = Q\_m + \mathbb{C}\_{gb} \left( \dot{\theta}\_p - \dot{\theta}\_m \right) + \mathbb{K}\_{gb} \left( \theta\_p - \theta\_m \right) + Q\_c\\ Q\_p = \mathbb{C}\_{Q\_p} \left( \dot{\theta}\_p, AR \right) \rho D\_p^5 \dot{\theta}\_p^2\\ AR = \mathcal{V}\_d / D\_p \theta\_p\\ Q\_c = Q\_{c\max} \sin(n\_h n\_d \theta\_m) \end{cases} \tag{1}$$

the motor, respectively, is the propeller-resistant torque, is a gust-induced disturbance torque, is the motor torque, , is the cogging torque, and is the maximum cogging torque, is the pole pairs number, is the harmonic index of the cogging disturbance, is the nondimensional torque coefficient of the propeller, is the propeller advance ratio, is the propeller diameter, is the air density, is the UAV forward speed, while and are the stiffness and the damping of the mechanical coupling joint. where *J<sup>p</sup>* and *θp*, *Jm*, and *θ<sup>m</sup>* are the inertia and the angular position of the propeller and the motor, respectively, *Q<sup>p</sup>* is the propeller-resistant torque, *Q<sup>d</sup>* is a gust-induced disturbance torque, *Q<sup>m</sup>* is the motor torque, *Qc*, is the cogging torque, and *Qcmax* is the maximum cogging torque, *n<sup>d</sup>* is the pole pairs number, *n<sup>h</sup>* is the harmonic index of the cogging disturbance, *CQ<sup>p</sup>* is the nondimensional torque coefficient of the propeller, *AR* is the propeller advance ratio, *D<sup>p</sup>* is the propeller diameter, *ρ* is the air density, V*<sup>a</sup>* is the UAV forward speed, while *Kgb* and *Cgb* are the stiffness and the damping of the mechanical coupling joint.

where and , , and are the inertia and the angular position of the propeller and

**Figure 4.** Full electric propulsion system: (**a**) mechanical scheme; (**b**) three-phase PMSM schematics with shorted turns in phase a (one pole pair and accessible neutral point). **Figure 4.** Full electric propulsion system: (**a**) mechanical scheme; (**b**) three-phase PMSM schematics with shorted turns in phase a (one pole pair and accessible neutral point).

### *2.3. Model of the PMSM with ITSC fault 2.3. Model of the PMSM with ITSC Fault*

The mathematical modelling of PMSM with ITSC faults that will be used for this work has been previously presented and applied in literature, demonstrating it to be very accurate [13,42–44]. The mathematical modelling of PMSM with ITSC faults that will be used for this work has been previously presented and applied in literature, demonstrating it to be very accurate [13,42–44].

The occurrence of a short circuit in a stator phase causes an asymmetry in the motor magnetic flux, due to an additional circuit with insulation resistance , in which the The occurrence of a short circuit in a stator phase causes an asymmetry in the motor magnetic flux, due to an additional circuit with insulation resistance *R<sup>f</sup>* , in which the shorted current *i <sup>f</sup>* flows, Figure 4b:

shorted current flows, Figure 4b: Once indicated, the fault extension along the winding with the parameter = / (defined as the ratio of the number of shorted turns to the total one), the winding affected by ITSC can be split into an undamaged and a damaged part, Figure 4b. Hence, the electrical equations can be written as: Once indicated, the fault extension along the winding with the parameter *µ* = *N<sup>f</sup>* /*N* (defined as the ratio of the number of shorted turns to the total one), the winding affected by ITSC can be split into an undamaged and a damaged part, Figure 4b. Hence, the electrical equations can be written as:

$$\mathbf{V}\_{\rm xn} = \mathbb{R}\mathbf{I} + \mathbb{L}\frac{d}{dt}\mathbf{I} + \mathbf{e}.\tag{2}$$

 = ℝ + <sup>+</sup> . (2) In Equation (2), = [ − , − , − , 0 ] is the applied voltages vector, = , , , is the currents vector, e is the BEMF vector, which can be expressed as in Equation (3) (in Equation (3) and in the following equations, "s" and "c" briefly indicates "sine" and "cosine" functions): In Equation (2), *Vxn* = [*V<sup>a</sup>* − *Vn*, *V<sup>b</sup>* − *Vn*, *V<sup>c</sup>* − *Vn*, 0 ] T is the applied voltages vector, *I* = h *ia*, *i<sup>b</sup>* , *ic*, *i f* iT is the currents vector, e is the BEMF vector, which can be expressed as in Equation (3) (in Equation (3) and in the following equations, "s" and "c" briefly indicates "sine" and "cosine" functions):

$$\mathbf{e} = \begin{bmatrix} (1 - \mu)e\_a \\ e\_b \\ e\_c \\ \mu e\_a \end{bmatrix} = -\lambda\_m \dot{\theta}\_e \begin{bmatrix} (1 - \mu)\mathbf{s}(\theta\_\varepsilon) \\ \mathbf{s}(\theta\_\varepsilon - 2/3\pi) \\ \mathbf{s}(\theta\_\varepsilon + 2/3\pi) \\ \mu \mathbf{s}(\theta\_\varepsilon) \end{bmatrix} \tag{3}$$

where = is the electrical angle of the motor and is the rotor magnet flux linkage. The resistance and inductance matrixes are: where *θ<sup>e</sup>* = *ndθ<sup>m</sup>* is the electrical angle of the motor and *λ<sup>m</sup>* is the rotor magnet flux linkage. The resistance and inductance matrixes are:

$$\mathbb{R} = \begin{bmatrix} (1-\mu)\mathbb{R} & 0 & 0 & \mathbb{R}\_f \\ 0 & \mathbb{R} & 0 & 0 \\ 0 & 0 & \mathbb{R} & 0 \\ \mu\mathbb{R} & 0 & 0 & -\left(\mu\mathbb{R} + \mathbb{R}\_f\right) \end{bmatrix}, \quad \mathbb{L} = \begin{bmatrix} (1-\mu)^2 L & 0 & 0 & 0 \\ 0 & L & 0 & 0 \\ 0 & 0 & L & 0 \\ \mu^2 L & 0 & 0 & -\mu^2 L \end{bmatrix}, \tag{4}$$

tion resistance, given by:

in which *R* is the phase resistance, *L* is the phase self-inductance, and *R<sup>f</sup>* is the insulation resistance, given by:

$$R\_f = k\_{R\_f} R (1 - \mu) \,\tag{5}$$

with *kR<sup>f</sup>* a factor depending on the insulation material.

By considering an inter-turn fault, the PMSM torque can be calculated as [44]:

$$Q\_m = \left( (1 - \mu)e\_a i\_a + e\_b i\_b + e\_c i\_c - \mu e\_a i\_f \right) / \dot{\theta}\_{m\nu} \tag{6}$$

Since the neutral point voltage *V<sup>n</sup>* is not null when a short circuit occurs, it is convenient to reformulate the electrical equations by extrapolating *Vn*. Hence, by applying Kirchhoff laws to the circuit in Figure 4b and by substituting in Equation (2), we have:

$$\mathbf{V} = \mathcal{R}\mathbf{I} + \mathcal{L}\frac{d}{dt}\mathbf{I} + \mathbf{E},\tag{7}$$

where the reformulated BEMF vector *E* is:

$$\mathbf{E} = \frac{1}{3} \begin{bmatrix} 2e\_{a\_f} - (e\_b + e\_c) \\ 2e\_b - \left(e\_{a\_f} + e\_c\right) \\ 2e\_c - \left(e\_{a\_f} + e\_b\right) \\ 3e\_f \end{bmatrix}'\tag{8}$$

The new resistance and inductance matrixes are:

$$
\mathcal{R} = \begin{bmatrix}
(1 - 2/3\mu)\mathcal{R} & 0 & 0 & 2/3\mathcal{R}\_f \\
\mu\mathcal{R}/3 & \mathcal{R} & 0 & -\mathcal{R}\_f/3 \\
\mu\mathcal{R}/3 & 0 & \mathcal{R} & -\mathcal{R}\_f/3 \\
\mu\mathcal{R} & 0 & 0 & -\left(\mu\mathcal{R} + \mathcal{R}\_f\right)
\end{bmatrix}, \quad \mathcal{L} = \begin{bmatrix}
(1 - 4/3\mu + 2/3\mu^2)L & 0 & 0 & 0 \\
(2 - \mu)\mu L/3 & L & 0 & 0 \\
(2 - \mu)\mu L/3 & 0 & L & 0 \\
\mu^2 L & 0 & 0 & -\mu^2 L
\end{bmatrix}, \tag{9}
$$

while *V* = [*Va*, *V<sup>b</sup>* , *Vc*, 0 ] T is a vector defining the voltage commands sent by the converter to the motor terminals (first three components) and the voltage on the shorted turns. Finally, since the PMSM is controlled via FOC technique, it is convenient to express the voltage vector via the Clarke–Park transformation, as a function of direct and quadrature voltages *V<sup>d</sup>* and *Vq*, generated by the closed-loop control algorithms:

$$
\begin{bmatrix} V\_a \\ V\_b \\ V\_c \\ V\_n \end{bmatrix} = \sqrt{\frac{2}{3}} \begin{bmatrix} \mathbf{c}(\theta\_\ell) & -\mathbf{s}(\theta\_\ell) \\ \mathbf{c}(\theta\_\ell - 2/3\pi) & -\mathbf{s}(\theta\_\ell - 2/3\pi) \\ \mathbf{c}(\theta\_\ell + 2/3\pi) & -\mathbf{s}(\theta\_\ell + 2/3\pi) \\ 0 & 0 \end{bmatrix} \begin{bmatrix} V\_d \\ V\_q \end{bmatrix} \tag{10}
$$

## **3. Fault-Tolerant Control System**

*3.1. FDI Algorithm Conceptualization*

When a PMSM operates at constant speed without faults, the current phasor in the Clarke plane draws a circular trajectory, while it follows an elliptical trajectory when an ITSC fault occurs [28–31]. Torque oscillations appear, at the motor level, with amplitudes that depend on the short-circuit extension.

This section aims to demonstrate that the geometrical parameters of this elliptical trajectory (major, minor axes length, and axes inclination) are univocally related to the location and extension of the ITSC. In particular, we will demonstrate that:

It is worth noting that, during constant speed operations, any elliptical trajectory of the motor current phasor in the Clarke plane (including the circular one, as a special case) can be reconstructed via the Fortescue decomposition [45], as the sum of three symmetrical and balanced rotating systems: *positive*-, *negative*-, and *zero-sequence* systems. *Aerospace* **2022**, *8*, x FOR PEER REVIEW 8 of 25

> When an ITSC fault occurs, the central point of the Y-connection is isolated, so the *zerosequence* current is zero, and the ellipse tracked by the current phasor (ˆ*I*) can be considered as the sum of two counter-rotating phasors: *positive*- (ˆ*I* <sup>+</sup>) and *negative-sequence* phasors ( ˆ*I* −), as given by: When an ITSC fault occurs, the central point of the Y-connection is isolated, so the *zero-sequence* current is zero, and the ellipse tracked by the current phasor () can be considered as the sum of two counter-rotating phasors: *positive-* () and *negative-sequence* phasors (), as given by:

$$\hat{I} = |\hat{I}(t)|e^{j\varphi(t)} = |\hat{I}^+|e^{j(\dot{\theta}\_\ell t + \varphi^+)} + |\hat{I}^-|e^{-j(\dot{\theta}\_\ell t - \varphi^-)}\,\,\,\,\,\,\tag{11}$$

where *ϕ*, *ϕ* <sup>+</sup>, and *ϕ* − are the phase angles of the resultant, *positive*-, and *negative*-*sequence* current phasors, respectively. As illustrated in Figure 5, the semi major and semi minor axis lengths are obtained by: Where , , and are the phase angles of the resultant, *positive-,* and *negative*-*sequence* current phasors, respectively. As illustrated in Figure 5, the semi major and semi minor axis lengths are obtained by:

$$\mathbf{s}\_{M,m} = |\mathbf{f}^+| \pm |\mathbf{f}^-|. \tag{12}$$

**Figure 5.** Decomposition of the current phasor in the Clarke plane into *positive-* and *negative-sequence* components. **Figure 5.** Decomposition of the current phasor in the Clarke plane into *positive*- and *negativesequence* components.

While the inclination of the ellipse major axis corresponds to the angle that maximizes the current phasor amplitude, up to: While the inclination of the ellipse major axis corresponds to the angle that maximizes the current phasor amplitude, up to:

$$\max\left(|\hat{I}|\right) = |\hat{I}^+| + |\hat{I}^-|.\tag{13}$$

The condition in Equation (13) is reached when the *positive-* and *negative-sequence* phasors have the same angular orientation: = ̇ |(||) + = −̇ |(||) + , (14) The condition in Equation (13) is reached when the *positive*- and *negative-sequence* phasors have the same angular orientation:

$$\left.\dot{\varphi}\_{\varepsilon l} = \dot{\theta}\_{\varepsilon} t\right|\_{\max(\left|\dot{I}\right|)} + \left.\varphi^{+} = -\dot{\theta}\_{\varepsilon} t\right|\_{\max(\left|\dot{I}\right|)} + \left.\varphi^{-}\right|\_{\max(\left|\dot{I}\right|)}\tag{14}$$

which leads to: which leads to:

m leads to:

$$
\varphi\_{el} = \left(\boldsymbol{\varrho}^{-} + \boldsymbol{\varrho}^{+}\right)/2.\tag{15}
$$

*positive-* and *negative-sequence* phasors, expressed by: <sup>∗</sup> = atan <sup>∗</sup> <sup>∗</sup> ⁄ , (16) To identify the ellipse inclination, it is thus crucial to evaluate the phase angles of the *positive*- and *negative-sequence* phasors, expressed by:

By introducing a complex analysis notation, the phasor ∗ can be also represented in

$$\boldsymbol{\sigma}^\* = \operatorname{atan} \left( i\_{\boldsymbol{\beta}}^\* / i\_{\boldsymbol{\alpha}}^\* \right) . \tag{16}$$

terms of a balanced three-phase space vector system, as:

the and axes in the Clarke plane.

where "∗" stands for "+" or "−", and *i* ∗ *<sup>α</sup>* and *i* ∗ *β* are the projections of the phasor ˆ*I* ∗ on the *α* and *β* axes in the Clarke plane.

By introducing a complex analysis notation, the phasor ˆ*I* ∗ can be also represented in terms of a balanced three-phase space vector system, as:

$$\hat{I}^\* = \sqrt{2/3} \left( \text{Re} \left( \hat{I}\_a^\* \right) + \text{Re} \left( \hat{I}\_b^\* \right) e^{j2\pi/3} + \text{Re} \left( \hat{I}\_c^\* \right) e^{j4\pi/3} \right) \tag{17}$$

Then, by expressing the space vector ˆ*I* ∗ into its real (*i* ∗ *α* ) and imaginary (*i* ∗ *β* ) components, we have:

$$\begin{split} \mathbf{i}\_a^\* &= \sqrt{2/3} \left[ \mathrm{Re} \left( \mathbf{\hat{I}}\_a^\* \right) - 1/2 \left( \mathrm{Re} \left( \mathbf{\hat{I}}\_b^\* \right) + \mathrm{Re} \left( \mathbf{\hat{I}}\_c^\* \right) \right) \right] = \sqrt{3/2} \mathrm{Re} \left( \mathbf{\hat{I}}\_a^\* \right), \\ \mathbf{i}\_\beta^\* &= \sqrt{2/3} \left[ \sqrt{3}/2 \left( \mathrm{Re} \left( \mathbf{\hat{I}}\_b^\* \right) - \mathrm{Re} \left( \mathbf{\hat{I}}\_c^\* \right) \right) \right] = \ast - \sqrt{3/2} \mathrm{Im} \left( \mathbf{\hat{I}}\_a^\* \right), \end{split} \tag{18}$$

The phase angles of the *positive*- and *negative-sequence* phasors with ITSC faults can be thus calculated via the *positive*- and *negative-sequence* phasors related to phase *a* only ( ˆ*I* + *<sup>a</sup>* and ˆ*I* − *a* ), given by the Fortescue transform:

$$
\begin{bmatrix}
\hat{I}\_a^+ \\
\hat{I}\_a^- \\
0 \\
\hat{I}\_f
\end{bmatrix} = \frac{1}{3} \begin{bmatrix}
1 & e^{j2\pi/3} & e^{j4\pi/3} & 0 \\
1 & e^{j4\pi/3} & e^{j2\pi/3} & 0 \\
1 & 1 & 1 & 0 \\
0 & 0 & 0 & 3
\end{bmatrix} \begin{bmatrix}
\hat{I}\_a \\
\hat{I}\_b \\
\hat{I}\_c \\
\hat{I}\_f
\end{bmatrix} \tag{19}
$$

The vector of phasors ˆ*I* = h ˆ*Ia*, ˆ*Ib* , ˆ*Ic*, ˆ*If* i can be obtained from:

$$\mathbf{J} = \left(\mathcal{R} + j\dot{\theta}\_{\varepsilon}\mathcal{L}\right)^{-1}(\mathcal{V} - \mathbf{E}),\tag{20}$$

where the voltage and BEMF phasors are:

$$
\hat{\mathbf{V}} = \sqrt{\frac{2}{3}} (V\_d + jV\_q) \begin{bmatrix} 1 \\ \begin{pmatrix} -1/2 - j\sqrt{3}/2 \\ -1/2 + j\sqrt{3}/2 \\ 0 \end{pmatrix} \\\tag{21}
$$

$$
\hat{\mathbf{E}} = j\lambda\_m \dot{\theta}\_e \begin{bmatrix} 1 - 2/3\mu \\ \mu - 1/2 - j\sqrt{3}/2 \\ \mu - 1/2 + j\sqrt{3}/2 \\ \mu \end{bmatrix} \tag{22}
$$

The solution of Equation (20) in terms of current phasors are presented by the polarcoordinate plot in Figure 6a, in which, given along the radius the ITSC fault extension (*µ* = *N<sup>f</sup>* /*N*), the inclination of the ellipse major axis of the current phasor trajectory is reported along the phase angle for each motor phase. In addition, Figure 6b shows the ratio of ellipse semi-axes lengths as a function of the ITSC fault extension parameter *µ*. It is worth noting that, if an ITSC fault occurs on phase *a* (*b*, *c*), the ellipse inclination will be 0 ◦ (240◦ , 120◦ ) ±15◦ .

**Figure 6.** Inclination of the ellipse major axis (**a**) and ellipses axis ratio (**b**) under ITSC faults of dif‐ ferent extensions and locations. **Figure 6.** Inclination of the ellipse major axis (**a**) and ellipses axis ratio (**b**) under ITSC faults of different extensions and locations.

The above‐mentioned discussion thus supports the development of an FDI algorithm based on the estimation of the geometrical characteristics of (generic) elliptical trajectories of current phasors in the Clarke plane, so that: The above-mentioned discussion thus supports the development of an FDI algorithm based on the estimation of the geometrical characteristics of (generic) elliptical trajectoriesof current phasors in the Clarke plane, so that:

The FDI concept essentially entails the identification of geometrical properties of an ellipse. Clearly, at a single monitoring time, a single point of the ellipse is known (ఈ, ఉ Figure 5); therefore, the fitting problem is under‐determined. On the other hand, the prob‐ lem becomes over‐determined if a sufficient number of measurements are used. The over‐ determined problem of fitting an ellipse to a set of data points arises in many applications such as computer graphics [46], hydraulic engineering [47], and statistics [48]. Many dif‐ ferent methods are proposed in literature for fitting the geometric primitives of an ellipse: voting/clustering methods (e.g., the Hough transform [49] and the RANSAC technique [50]) are robust but they have poor accuracy and require large computing resources, opti‐ mization methods and direct least‐square methods [51] are accurate, but, in case of non‐ linear applications such as an ellipse fitting, they apply iterative techniques. The FDI concept essentially entails the identification of geometrical properties of an ellipse. Clearly, at a single monitoring time, a single point of the ellipse is known (*iα*, *i<sup>β</sup>* Figure 5); therefore, the fitting problem is under-determined. On the other hand, the problem becomes over-determined if a sufficient number of measurements are used. The over-determined problem of fitting an ellipse to a set of data points arises in many applications such as computer graphics [46], hydraulic engineering [47], and statistics [48]. Many different methods are proposed in literature for fitting the geometric primitives of an ellipse: voting/clustering methods (e.g., the Hough transform [49] and the RANSAC technique [50]) are robust but they have poor accuracy and require large computing resources, optimization methods and direct least-square methods [51] are accurate, but, in case of nonlinear applications such as an ellipse fitting, they apply iterative techniques.

Here we instead use a direct least‐square method, as proposed by Halìř and Flusser [52], who, starting from the quadratically‐constrained least‐squares minimization postu‐ lated by Fitzgibbon [53], reformulated the problem via partitioning techniques, by obtain‐ ing a very robust and computationally‐effective algorithm. Here we instead use a direct least-square method, as proposed by Halìˇr and Flusser [52], who, starting from the quadratically-constrained least-squares minimization postulated by Fitzgibbon [53], reformulated the problem via partitioning techniques, by obtaining a very robust and computationally-effective algorithm.

Once given the vectorial conic definition of an ellipse, Once given the vectorial conic definition of an ellipse,

$$
\Gamma \cdot \gamma = 0,\tag{23}
$$

in which ൌ ሾଶ, , ଶ, , , 1ሿ and ൌ ሾ, , , , , ሿ are vectors defining the Carte‐ sian coordinates of the ellipse points and the ellipse coefficients, respectively, the over‐ determined fitting problem to a set of coordinate points and (where ൌ 1, … , , in which **Γ** = - *α* 2 , *αβ*, *β* 2 , *α*, *β*, 1 and *γ* = [*A*, *B*, *C*, *D*, *E*, *F*] T are vectors defining the Cartesian coordinates of the ellipse points and the ellipse coefficients, respectively, the overdetermined fitting problem to a set of coordinate points *α<sup>i</sup>* and *β<sup>i</sup>* (where *i* = 1, . . . , *n*, and

and is greater than the number of ellipse coefficients, i.e., >6) can be solved by the

following eigenvalues problem (for more details, see Appendix A):

*n* is greater than the number of ellipse coefficients, i.e., *n* > 6) can be solved by the following eigenvalues problem (for more details, see Appendix A):

$$\begin{cases} \mathbb{M}\gamma\_1 = \lambda\gamma\_1\\ \begin{array}{l} \gamma\_1^\mathrm{T} \mathbb{C}\_1 \gamma\_1 = 1\\ \gamma\_2 = -\mathbb{S}\_3^{-1} \mathbb{S}\_2^\mathrm{T} \gamma\_1\\ \gamma = \left(\gamma\_1 \ \gamma\_2\right)^\mathrm{T} \end{array} \tag{24}$$

where M is the reduced scatter matrix,

$$\mathbb{M} = \mathbb{C}\_1^{-1} \left( \mathbb{S}\_1 - \mathbb{S}\_2 \mathbb{S}\_3^{-1} \mathbb{S}\_2^{\mathrm{T}} \right), \tag{25}$$

$$\mathbb{S}\_1 = \mathbb{D}\_1^\mathrm{T} \mathbb{D}\_1, \quad \mathbb{S}\_2 = \mathbb{D}\_1^\mathrm{T} \mathbb{D}\_2, \quad \mathbb{S}\_3 = \mathbb{D}\_2^\mathrm{T} \mathbb{D}\_2,\tag{26}$$

$$\mathbb{D}\_1 = \begin{bmatrix} \alpha\_1^2 & \alpha\_1 \beta\_1 & \beta\_1^2 \\ \vdots & \vdots & \vdots \\ \alpha\_n^2 & \alpha\_n \beta\_n & \beta\_n^2 \end{bmatrix}, \begin{bmatrix} \mathbb{D}\_2 = \begin{bmatrix} \alpha\_1 & \beta\_1 & 1 \\ \vdots & \vdots & \vdots \\ \alpha\_n & \beta\_n & 1 \end{bmatrix} \end{bmatrix} \tag{27}$$

$$\mathbb{C}\_1 = \begin{bmatrix} 0 & 0 & 2 \\ 0 & -1 & 0 \\ 2 & 0 & 0 \end{bmatrix} \, \tag{28}$$

$$\boldsymbol{\gamma} = \begin{bmatrix} \boldsymbol{\gamma} \ \boldsymbol{\gamma} \end{bmatrix}^{\mathrm{T}} \tag{29}$$

with the ellipse coefficient vector segmented in *γ*<sup>1</sup> = - *A B C*<sup>T</sup> and *γ*<sup>2</sup> = - *D E F*<sup>T</sup> . The solution of Equation (24) corresponds to the eigenvector *γ* yielding a minimal nonnegative eigenvalue *λ*. Once obtained *γ*, the lengths of major and minor semi-axes *s<sup>M</sup>* and *s<sup>m</sup>* are given by [54]:

$$s\_{M,m} = \frac{\sqrt{\frac{2(AE^2 + \mathbb{C}D^2 - BDE + (B^2 - 4AC)F)}{\left(A + \mathbb{C} \pm \sqrt{(A - \mathbb{C})^2 + B^2}\right)^{-1}}}}{4AC - B^2}}{4AC - B^2} \tag{30}$$

while the major axis inclination *ϕel* is [54]:

$$\varphi\_{el} = \begin{cases} 0 & \text{for } \mathcal{B} = 0, \ A < \mathcal{C} \\ \pi/2 & \text{for } \mathcal{B} = 0, \ A > \mathcal{C} \\ 1/2 \cot^{-1}((A - \mathcal{C}) \text{ / } \mathcal{B}) & \text{for } \mathcal{B} \neq 0, \ A < \mathcal{C} \\ \pi/2 + 1/2 \cot^{-1}((A - \mathcal{C}) \text{ / } \mathcal{B}) & \text{for } \mathcal{B} \neq 0, \ A > \mathcal{C} \end{cases} \tag{31}$$

## *3.2. FDI Algorithm Design and Implementation*

One of the basic limitations of the proposed FDI concept is that the technique is defined by assuming constant speed motor operations. During unsteady conditions, the current phasor changes its amplitude while rotating in the Clarke plane. Consequently, the ellipse identification does not work appropriately, and it could cause false alarms.

To overcome the problem, it is worth noting that, differently from the constant speed operations with ITSC faults, the axes of the ellipse reconstructed in unsteady conditions rotate. For this reason, the developed FDI algorithm detects an ITSC fault only if the value of the ellipse major axis is constant and stationary. This FDI logic has been implemented as represented by the flow chart in Figure 7: at each *k*-th monitoring sample, if the Boolean variable *mon*(*k*) is true (Equation (32)), a fault counter *n* (*k*) *c* is increased by two, otherwise it is reduced by one. When the fault counter reaches a predefined maximum value (*nth*), the algorithm outputs a true Boolean fault flag *fx*, which indicates that a fault on the phase *x* (= *a*, *b*, *c*) is detected and isolated.

$$
gamma^{(k)} = \left(\Delta\_d^{(k)} \ge \varepsilon\_{\rm{dtlh}} \land \ \Delta\_i^{(k)} \le \varepsilon\_{\rm{itlh}}\right) \tag{32}
$$

where:

$$\Delta\_d^{(k)} = \left| s\_M^{(k)} - s\_m^{(k)} \right|,\tag{33}$$

$$
\Delta\_i^{(k)} = \min \left( \left. \Delta\_i^{(k)} \right|\_{a'} \left. \Delta\_i^{(k)} \right|\_{c'} \left. \Delta\_i^{(k)} \right|\_{b} \right), \tag{34}
$$

$$
\left| \Delta\_i^{(k)} \right|\_a = \left| \varphi\_{el}^{(k)} \right|\_\prime \ \Delta\_i^{(k)} \Big|\_b = \left| \varphi\_{el}^{(k)} \right| - \frac{2\pi}{3} \ \ \Delta\_i^{(k)} \Big|\_c = \left| \varphi\_{el}^{(k)} \right| + \frac{2\pi}{3} \ \tag{35}
$$

**Figure 7.** Flow chart for the FDI logic. **Figure 7.** Flow chart for the FDI logic.

where: Δ () <sup>=</sup> () <sup>−</sup> () , (33) The details and a dedicated discussion about the FDI parameters design are proposed in Section 4, while tables reporting the FDI parameters are reported in Appendix B.

## *3.3. Fault Accommodation Algorithm*

Δ () <sup>=</sup> min Δ () , Δ () , Δ () , (34) Δ () = () , Δ () = () <sup>−</sup> 2 <sup>3</sup> , Δ () = () + 2 <sup>3</sup> , (35) The details and a dedicated discussion about the FDI parameters design are proposed in Section 4, while tables reporting the FDI parameters are reported in Appendix B. After the FDI algorithm detects and isolates the ITSC fault, the motor still has two active phases. To maintain the performances after the fault, the FTC should accommodate the system by restoring the operation of the current phasor in the Clarke plane as it was in healthy conditions. This is here achieved by enabling the control of the central point of the motor Y-connection, via the stand-by leg of the four-leg converter, and by applying an RCFTC technique [9,55].

*3.3. Fault Accommodation Algorithm* After the FDI algorithm detects and isolates the ITSC fault, the motor still has two active phases. To maintain the performances after the fault, the FTC should accommodate the system by restoring the operation of the current phasor in the Clarke plane as it was To maintain the torque performances with only two phases, the currents amplitudes in the healthy phases in the rotor reference frame (*i* # *x f* , *i* # *y f*) must increase by <sup>√</sup> 3 w.r.t. those before the fault, and they must shift by 60◦ along the electrical angle. In addition, the amplitude of the current flowing into the central point (*i* # *n* ) must be <sup>√</sup> 3 times those in the healthy phases:

$$\begin{cases} i\_{xf}^{\#} = \sqrt{2} \left( i\_d^{\#} \text{c} \left[ \theta\_\varepsilon + \frac{2\pi}{3} \left( m + \frac{7}{4} \right) \right] - i\_q^{\#} \text{s} \left[ \theta\_\varepsilon + \frac{2\pi}{3} \left( m + \frac{7}{4} \right) \right] \right) \\ i\_{yf}^{\#} = \sqrt{2} \left( i\_d^{\#} \text{c} \left[ \theta\_\varepsilon + \frac{2\pi}{3} \left( m + \frac{5}{4} \right) \right] - i\_q^{\#} \text{s} \left[ \theta\_\varepsilon + \frac{2\pi}{3} \left( m + \frac{5}{4} \right) \right] \right) \\\ i\_{wf}^{\#} = 0 \\\ i\_n^{\#} = \sqrt{6} \left( i\_d^{\#} \text{c} \left( \theta\_\varepsilon + \frac{2\pi}{3} m \right) - i\_q^{\#} \text{s} \left( \theta\_\varepsilon + \frac{2\pi}{3} \right) \right) \end{cases} \tag{36}$$


*2.4. Numerical Simulation*

The numerical analysis software COMSOL multiphysics5.5 was used for simulation analysis in this study. The laminar flow module and particle trajectory tracking module in the software were used to simulate the separation process. The particle trajectory tracking module couples the effect of multiple force fields on the particle trajectory, calculating the trajectory of particles in the channel.

First, a 2D sketch of the microfluidic chip was drawn by AutoCAD2021 and then imported into COMSOL multiphysics5.5. The model needed to be meshed before the analysis, and the meshing was finer to achieve a better simulation. According to the actual situation of the subject, the Reynolds number did not exceed 1000, so the fluid was set as laminar flow. The Reynolds numbers in collection area 1 and 2 were 173 and 346, respectively. The particle density was 1.05 g cm−3 and the aerodynamic diameter was 0.5~8 μm. The inlet flow rate was 12.5 mL min−1, and the wall was set to be slipproof. During the simulation, the collisions between the particles and the walls of the microfluidic device were inelastic, and the particles stuck to the walls to calculate the particle collection rate. Particle tracking uses Newton's law of motion to solve differential equations. The research object of this topic is the movement of spores in the air. The spores were mainly driven by drag force in the microfluidic chip. The particle tracking module was used to simulate the particle trajectory. The drag force satisfies Stokes' law, and the temperature and absolute pressure are standard states. Then, the simulation boundary was, the corresponding entrance and exit were set, and 100 particles were released at the entrance each time to ensure the reliability of particle collection efficiency. Since the simulation of this study needs to study the laminar flow and particle motion state in the microfluidic chip, transient and steady-state solvers were configured. The steady-state solver was used to study laminar flow, and the transient solver was used to study the particle. In the case of motion, the steady-state calculation amount is small, and one can choose to solve it directly. The solver can choose PARDISO to meet the requirements, and then choose the multi-threaded nested analysis pre-sorting algorithm, and the configuration of the steady-state solver is completed. The transient solver is more complicated to solve, and it needs to be solved iteratively. The GMRES solver was selected to solve the trajectory of the particle.