Finite Element Analysis of Laminar Heat Transfer within an Axial-Flux Permanent Magnet Machine

Abstract

Axial-Flux Permanent Magnet (AFPM) machines have gained popularity over the past few years due to their compact design. Their application can be found, for example, in the automotive and medical sectors. For typically considered materials, excessive heat can be generated, causing possible irreversible damage to the magnets, bonding, or other structural parts. In order to optimize cooling, knowledge of the flow and the consequent temperature distribution is required. This paper discusses the flow types and heat transfer present inside a typical AFPM machine. An Isogeometric Analysis (IGA) laminar-energy model is developed using the Nutils open-source Python package. The developed analysis tool is used to study the effects of various important design parameters, such as the air-inlet, the gap-length, and the rotation speed on the heat transfer in an AFPM machine. It is observed that the convective heat transfer at the stator core is negatively affected by adding an air-inlet. However, the heat dissipation of the entire stator improves as convective heat transfer occurs within the air-inlet.

1. Introduction

Rotating disks are a common geometry in a variety of engineering applications. Well-known examples include radial compressors, turbines, disk brakes, friction pumps, and electrical machinery. This manuscript focuses on electrical machinery. More specifically, the electrically powered Axial-Flux Permanent Magnet (AFPM) machine is the main application of interest. In general, the AFPM machine consists of two axially centered disks referred to as the rotor and stator. The rotor is composed of a permanent magnet array that is attached to a back-iron. The stator is composed of current-carrying windings or coils and a soft-magnetic core. Torque is produced through the interaction of the magnetic field created by the permanent magnets and the magnetic field generated by the windings. As a result of a commutation strategy, the rotor is put into motion. A typical AFPM machine is shown in Figure 1, which has two rotors on either side of the stator. The current-carrying coils or windings are located on the stator (not visible in the figure). The back-iron plate is in reality mounted to a specific object that is required to be put into motion, e.g., a wheel, which typically encloses the majority of the machine.

The interest in brushless permanent magnet (PM) machines, of which the AFPM is an example, has been increasing over the past two decades [1]. The main reason for this is the higher efficiency and power density compared to induction and brushed motors. AFPM machines are more compact than their common Radial-Flux PM (RFPM) counterparts. Applications of AFPM machines range from the high-performance automotive sector to the medical sector and even consumer electronics. This diversity in applications creates additional design problems.

A design problem specifically related to the AFPM machine is the iron stator material cost. For small high-performance motors, expensive soft-magnetic composite materials (SMC) are used that exhibit negligible eddy current losses. Medical equipment, such as CT-scanners, require large material volumes, however, and therefore favor the use of cheaper electrical steel. Due to the higher electrical conductivity in cheap iron materials, a substantial amount of eddy currents is induced by the alternating magnetic field due to the rotation of the magnets. These eddy currents are dissipated in the form of heat. Because of the compact design of AFPM machines, this generated heat becomes challenging to control. Consequently, temperature-dependent material properties such as electrical and thermal conductivity or even the remanent magnetic flux density of the permanent magnets are affected, decreasing the motor’s performance.

To aid in the design of AFPM machines, an electromagnetic finite element model was recently proposed by Friedrich et al. [2] to efficiently analyze the eddy currents reduction in the stator core. Adjustments in material properties and geometric design, i.e., adding slits in the circumferential and radial directions, are considered. To analyze the above-mentioned temperature dependent effects in AFPM machines, the model of Ref. [2] needs to be enriched with thermal modeling capabilities. On the one hand, this enrichment pertains to the inclusion of the fluid flow through the AFPM machine, as this flow plays a critical role in the cooling of the machine. On the other hand, a thermal energy balance must be introduced to describe the temperature field in the AFPM machine. The state of the art of both these aspects in the context of this manuscript will be reviewed in Section 1.1 and Section 1.2, respectively. The research objective is subsequently introduced in Section 1.3.

1.1. Flow Characteristics

Research on rotating disk flows, as encountered in the AFPM machine application, increased drastically after the pioneering Ph.D. research of Ekman in 1905 [3], who investigated the air and sea flow behavior due to earth’s rotation. Although this research was conducted in the field of oceanography, its application can be extended to geophysics, astrophysics, and turbomachinery. Apart from localized differences, it has been shown that the general flow field in an AFPM machine is similar to a simplified rotor-stator configuration [4,5,6] to which the generic theory of Ref. [3] applies. A schematic representation of such a rotor-stator configuration simplification is visualized in Figure 2. In this configuration, two coaxial parallel disks, positioned at a specific distance, rotate with respect to each other with a rotation speed $\Omega$. The disks are separated by an axial gap-length d. The gap may be enclosed by an additional static external cylinder (shroud), which alters the flow behavior as will be discussed later. The imposed difference -stator radii caused by the shroud is often negligible and hence these radii are typically both set to the same radius R. A centrifugal throughflow Q can be superposed through an axially centered inlet with radius $R_i$.

A variant of this problem was first investigated by only considering the rotor (excluding the stator parts) by von Kármán [7]. He showed that with an infinite radius and laminar flow, the flow is confined in a thin boundary layer on the disk, moving outwards due to centrifugal forces. Consequently, an axial inflow is obtained, satisfying the continuity equation. He also studied the turbulent flow-case by means of momentum integral methods. The von Kármán analysis was also followed by Bödewadt [8], who numerically investigated a rigid-body rotating flow over a stationary plane, i.e., a stator. Both these phenomena, i.e., the outward moving boundary layer and the rigid-body rotation, were observed for the flow between two infinite parallel disks, as considered by Batchelor who solved the differential equations relative to the stationary axisymmetric flow [9]. He noticed the formation of a nonviscous core in the solid body rotation, confined between the two boundary layers which develop on the disks. These boundary layers on both the rotor and the stator are, respectively, referred to as Ekman or von Kármán and Bödewadt layers.

Figure 3a–c illustrates the Batchelor flow, where the radial velocity component $u_r$ in Figure 3b clearly shows three distinct zones. A rotor boundary layer is established due to the centrifugally (outward) pumped center fluid. Consequently, an opposite radially inward flow is found at the stator, moving to the center. The core-fluid itself does not move in the radial direction but merely in the tangential or angular direction, as shown in Figure 3a. The ratio between the radial and tangential core velocity is known as the core-swirl ratio and is quantified by an entrainment coefficient [10]. The distinction of the three zones by Batchelor as depicted in Figure 3a–c, was a subject of debate from 1953 to 1983 between Batchelor and Stewartson. Stewartson performed similar analyses but did not observe the core swirl observed by Batchelor [11]. The absence of a core swirl was further confirmed by numerical solutions of the Navier–Stokes equations. Stewartson reasoned that there was no need for a boundary layer on the stator and that the tangential velocity component is equal to zero everywhere, except near the rotor, as illustrated in Figure 3d. An outward-pointing radial velocity is then obtained, as exemplified in Figure 3e.

The debate between Batchelor and Stewartson was settled by Zandbergen and Dijkstra [13], who showed that the solutions derived by means of the similarity method are not unique. This means that both solutions provided by Batchelor and Stewardson are possible, depending on the flow conditions.

To describe these flow conditions, three dimensionless parameters are defined, viz., the aspect ratio G, the throughflow coefficient $C_w$ and the rotational or global Reynolds number Re:

$$\mathrm{G}=\frac{d}{R},C_w=\frac{Q}{\nu R},\mathrm{and}\mathrm{Re}=\frac{\Omega R^2}{\nu }.$$

(1)

In these dimensionless numbers, d [m] is the axial gap and R [m] the outer radius of the rotor (see Figure 2). Q [m${}^3$/s] represents the superimposed throughflow, $\Omega$ [rad/s] the angular velocity of the rotor disk and $\nu$ [m${}^2$/s] the kinematic viscosity. It is common practice to express the Reynolds number using either the local radial position r or the axial gap d as the characteristic length. In the remainder, such use will be explicitly denoted by either a subscript r or d. On account of the above-defined aspect ratio and global Reynolds number it follows that ${\mathrm{Re}}_d=\mathrm{Re}{\mathrm{G}}^2$. It is noted that whether or not a shroud is present, respectively, referred to as either a closed- or open-end configuration, also influences the flow behavior.

Brady and Durfelofsky [14] investigated the laminar flow of both closed- and open-end rotor-stator configurations with a large but finite radius. By means of an asymptotic-numerical analysis they noticed that Batchelor and Stewartson flows are, respectively, found in closed- and open-end rotor-stator configuration for the case of ${\mathrm{Re}}_d>80$. At lower Reynolds numbers, ${\mathrm{Re}}_d\le 80$, both open- and closed-end flows show no significant difference as they correspond to the Batchelor flow type near the axis of rotation. A numerical and experimental study on the turbulent flow case [15] showed that either a Batchelor or Stewartson flow is obtained depending on the superposed throughflow in a closed-end configuration for ${10}^6\le \mathrm{Re}\le 4.15\times {10}^6$. When dealing with no or a centripetal (radially inward) throughflow, a Batchelor flow was observed. In contrast, a centrifugal (radially outward) imposed throughflow resulted in a Stewartson flow. This observation agrees with the positive radial velocity observed in Figure 3e, which is induced by the positive radial outflow. It should, however, be noted that a relatively low centrifugal throughflow creates a case in which both Stewartson and Batchelor flows can coexist, which stipulates that there exist cases for which no clear distinction between both flow types can be made.

Rotor-stator fluid structures have been numerically analyzed extensively, with the main focus on turbulent flows [15,16,17]. A complicating factor in these analyses is the presence of a transition region from laminar to turbulent flow [18,19]. Moreover, under certain flow conditions a three-dimensional (3D) transient S-vortex can occur [20]. These additional complexities in the flow are restricted to the turbulent flow cases.

1.2. Thermal Characteristics

To model the thermal behavior of AFPM machines, besides the cooling flow also the heat transfer must be taken into account. A heat source emanates from the Joules and eddy current losses in the coils, and from the eddy current, excess and hysteresis losses in the core [2,21]. In the remainder, only eddy currents are considered as the source of heat production. The resulting temperature field and heat transfer are the primary interests of this work as these affect the overall thermal performance of the AFPM machine. More specifically, they influence the remanent magnetization of the magnet and the electrical conductivity of the materials. Additionally, it is important to ensure that the machine remains below a specific maximum temperature. Surpassing such a temperature limit will affect the structural integrity and could potentially demagnetize the permanent magnets.

The heat transfer is quantified using the Nusselt number, which in general is defined as

$${\mathrm{Nu}}_{\mathrm{R}}=\frac{h_{c}R}{k},$$

(2)

in which $h_c$ [W/m${}^2$K] is the convective heat transfer coefficient, R [m] the outer radius, and k [W/mK] the conduction coefficient of the fluid. The Nusselt number relates the convective heat transfer caused by the fluid motion to the conductive heat transfer. From the perspective of the thermal design, higher values for the Nusselt number are desired, as these correspond to an improved heat transfer, and hence improved cooling characteristics. Concerning the AFPM machine, the Nusselt number can be used to evaluate the effect of changing the axial gap-length, the addition of cooling vanes, or even forced throughflow.

Although many experiments and analyses have been conducted on the heat transfer along the rotor [10,22,23], research on the heat transfer at the stator side is limited. Yuan et al. [24] conducted both experiments and numerical analyses on an open-end rotor-stator system without throughflow. In their analyses they have considered both disks to be adiabatic, apart from the stator-gap surface, on which an isothermal heat distribution was applied. For a rotational Reynolds number range of $1.42\times {10}^5\le \mathrm{Re}\le 3.33\times {10}^5$, they concluded that there exists a gap-length which optimizes the average Nusselt number. The optimal gap-length was found to decrease with increasing Reynolds number. Related to this optimum gap-length is the occurrence of a recirculating flow zone near the periphery of the system. Although only visible at relatively large gap-lengths, this circulation zone results in a decrease of the Nusselt number. This effect has also been observed in Ref. [22], who used a Direct Numerical Simulation (DNS) to model the heat transfer within two counter-rotating disks. In contrast, when such a recirculation is not present, a steep increase in the local Nusselt number near the periphery is observed. This observation is in agreement with the study presented in Ref. [25]. For the range of $0.0106\le \mathrm{G}\le 0.0297$ and $3.7\times {10}^4\le \mathrm{Re}\le 5.6\times {10}^5$, the observed increase in local Nusselt number was attributed to the ingress of ambient fluid along the stator as a result of the rotor’s pumping effect (see Figure 3b).

1.3. Research Objective

Despite the fact that the thermal characteristics of AFPM machines are a critical design aspect, the current understanding of these characteristics is limited. Further research is particularly needed with respect to the heat transfer at the stationary disk of open-end rotor-stator configurations. In this manuscript we present a computational analysis of the thermal behavior of rotor-stator configurations that are representative for a broad class of AFPM machines. Compared to the results available in the literature, our analysis improves the insight into the thermal characteristics of AFPM machines on account of two fundamental model enhancements. First, the developed model allows the stator-gap surface to be heated non-uniformly, which is an essential model feature required to represent the spatial distribution of the eddy current losses in AFPM machines. Second, the model is capable of solving the temperature distribution inside both the solid and the fluid regions, which enables the representation of heat dissipation across all interfaces. To discretize our model, the Isogeometric Analysis (IGA) paradigm is considered [26]. This finite-element-type discretization technique employs spline basis functions to represent both the geometry of the model and the physical fields. IGA has been demonstrated to be a suitable approximation technique for a wide range of problems, including problems in fluid flow analysis (e.g., [27,28,29]) and thermal analysis (e.g., [30,31]). In our work we particularly leverage the advantageous properties of IGA with respect to approximation properties for mixed formulations [32,33,34] and (local) refinement capabilities [35,36,37,38]. To focus on the thermal performance of AFPM machines, we restrict ourselves to the laminar flow case. Our simulations give additional insights into the effect of design improvements on the thermal performance of AFPM machines. The extent to which these insights can be transferred to the turbulent flow case is discussed.

The structure of this paper is as follows: in Section 2 the finite element flow model is introduced and validated using several benchmark problems. In Section 3, the thermal model is considered and validated in terms of temperature and heat transfer quantification, i.e., the Nusselt number. The flow model and thermal model are then combined to investigate the AFPM machine characteristics in Section 4, while detailed analysis results are presented in Section 5. Finally, conclusions and recommendations are presented in Section 6.

2. Flow Model

The flow model is developed in the framework of the finite element method (FEM), making use of the Python-based open-source FEM library Nutils [39]. The Isogeometric analysis (IGA) paradigm introduced by Hughes and co-workers [26] is employed. The control of inter-element regularity in IGA enables it to produce globally continuous and smooth fields. For a wide range of cases, in particular, for problems with smooth solutions, IGA has been demonstrated to yield a similar accuracy as traditional FEM using significantly fewer degrees of freedom (dofs) [40]. This benefit in terms of accuracy per dof generally translates into an improvement in terms of computational effort when IGA is considered instead of FEM. A detailed comparison of the required computational effort is tedious, however, as it also requires consideration of matrix assembly (incl., numerical integration) efforts, (iterative) solver performance, post-processing costs, differences in implementation, etc. Such a detailed comparison is considered beyond the scope of this manuscript.

In our work, special attention is given to the 2D axisymmetry assumption, which states that the gradient in the tangential direction is zero for any considered field. This assumption is often used-stator numerical analyses, due to favorable computational effort compared to three-dimensional simulations and its good agreement with experimental data in the turbulent regime [15,19,24]. A representative axisymmetric domain with the corresponding boundary conditions is visualized in Figure 4. Note that the considered model (formulation and discretization) can be extended ex mutatis mutandis to three dimensions, which poses no additional complexities concerning the isogeometric approach proposed in this paper. However, the treatment of a three-dimensional geometry will become demanding in terms of degrees of freedom and computational effort. The extension to three dimensions would require the efficient use of adaptive refinement techniques (for example, Refs. [41,42,43]), consideration of iterative solvers [44], and, most likely, the consideration of highly-scalable parallel computing libraries.

2.1. Strong Form

We consider an arbitrary closed 2D axisymmetric domain $\overline{\Omega }$, consisting of an open domain $\tilde{\Omega }$ bounded by a piecewise smooth boundary $\Gamma =\partial \tilde{\Omega }$ (see Figure 4). The viscous incompressible Navier–Stokes equations and mass conservation equation in steady-state on such a closed domain are given by

$$\begin{array}{cc}\hfill & \rho \mathbf{u}·\left(\nabla \mathbf{u}\right)=\nabla ·\left[2\mu {\nabla }^s\mathbf{u}-p\mathbf{I}\right]+\mathbf{f},\hfill \end{array}$$

(3a)

$$\begin{array}{cc}\hfill & \nabla ·\mathbf{u}=0,\hfill \end{array}$$

(3b)

with

$$\mathbf{u}=\mathbf{g}\mathrm{at}{\Gamma }_{\mathrm{D}}\mathrm{and}\mathbf{n}·\nabla \mathbf{u}=\mathbf{0}\mathrm{at}{\Gamma }_{\mathrm{N}}.$$

(4)

In the above equations, $\mathbf{u}$ [m/s] is the velocity vector, p [Pa] the pressure, $\mathbf{f}$ [N/m${}^3$] the body force, $\rho$ [kg/m${}^3$] the fluid density and $\mu$ [m${}^2$/s] the dynamic viscosity. A fixed velocity vector $\mathbf{g}$ [m/s] is assigned to the Dirichlet boundary. Note that it is also assumed that both the density and dynamic viscosity have a fixed value, i.e., there is no dependency on the temperature or pressure. The validity of this assumption has been investigated and found to be negligible (see Appendix B). The entire boundary is subdivided into a Dirichlet ${\Gamma }_{\mathrm{D}}$ and Neumann ${\Gamma }_{\mathrm{N}}$ part. The former is used as a wall boundary ${\Gamma }_w$, which can either represent a no-slip or rotating wall. In Equation (3a) ${\nabla }^s$ is the symmetric gradient, defined as

$${\nabla }^s\mathbf{u}=\frac{1}{2}\left(\nabla \mathbf{u}+\nabla {\mathbf{u}}^T\right).$$

(5)

Note the factor 2 in Equation (3a) is a consequence of this definition. It is important to note that the gradients in Equation (3) are taken with respect to the Cartesian reference frame. In our implementation, axisymmetry is modeled by mapping this coordinate system onto a cylindrical reference frame and then equating all circumferential derivatives to zero. The primary advantage of this implementation is that the flow model maintains its general form and could, as a result, be applied to non-axisymmetric flow problems without alterations.

2.2. Mixed Variational Form

In order to derive the weak formulation, appropriate spaces are first to be defined. Let $\left(u,p\right)\in \{\mathcal{S}\times \mathcal{P}\}$ and $\left(v,q\right)\in \{\mathcal{V}\times \mathcal{Q}\}$ denote the trial and test function spaces for the mixed problem, which are chosen to be identical

$$\begin{array}{cc}\hfill & \mathcal{S}:=\{\mathbf{u}\in H^1\left(\tilde{\Omega }\right)|\mathbf{u}·\mathbf{n}=0\mathrm{on}{\Gamma }_{rot}\},\hfill \end{array}$$

(6a)

$$\begin{array}{cc}\hfill & \mathcal{P}:=L^2\left(\tilde{\Omega }\right),\hfill \end{array}$$

(6b)

and

$$\begin{array}{cc}\hfill & \mathcal{V}:=\{\mathbf{v}\in H^1\left(\tilde{\Omega }\right)|\mathbf{v}·\mathbf{n}=0\mathrm{on}{\Gamma }_{rot}\},\hfill \end{array}$$

(7a)

$$\begin{array}{cc}\hfill & \mathcal{Q}:=L^2\left(\tilde{\Omega }\right),\hfill \end{array}$$

(7b)

with a no-penetration condition at the symmetry or axis of rotation boundary, ${\Gamma }_{rot}$, and with $L^2$ and $H^1$ denoting the space of square-integrable functions and the first order Sobolev space, respectively. Finite-dimensional versions of these spaces are defined as $\{{\mathcal{S}}^h\times {\mathcal{P}}^h\}\subset \{\mathcal{S}\times \mathcal{P}\}$ and $\{{\mathcal{V}}^h\times {\mathcal{Q}}^h\}\subset \{\mathcal{V}\times \mathcal{Q}\}$, where the superscript h is used to indicate that a space is finite dimensional.

For the Navier–Stokes equations, the required Sobolev space is of the first order. The velocity normal is strongly imposed on the rotation axis boundary ${\Gamma }_{rot}$ through the vector spaces defined in Equations (6) and (7). The wall boundary condition (stationary or rotating) is achieved by imposing the velocity constraints weakly by means of Nitsche’s method (see Section 2.2.1). The residual of the mixed Galerkin discretized variational formulation for an arbitrary domain $\tilde{\Omega }$, with a piecewise smooth boundary $\Gamma =\partial \tilde{\Omega }$ then states:

Find $\{{\mathbf{u}}^h,p^h\}\in {\mathcal{S}}^h\times {\mathcal{P}}^h$, such that:

$$\begin{array}{cc}\hfill & {\left({\mathbf{v}}^h,{\mathbf{u}}^h·\left[\nabla {\mathbf{u}}^h\right]\right)}_{\tilde{\Omega }}-{\left(\nabla ·{\mathbf{v}}^h,\frac{1}{\rho }{p}^h\right)}_{\tilde{\Omega }}+{\left({\nabla }^s{\mathbf{v}}^h,2\nu {\nabla }^s{\mathbf{u}}^h\right)}_{\tilde{\Omega }}\hfill \\ & -{\left({\mathbf{v}}^h,\frac{1}{\rho }{\mathbf{f}}^h\right)}_{\tilde{\Omega }}+{\left({\mathbf{v}}^h·\mathbf{n},p^h\right)}_{\Gamma }-{\left({\mathbf{v}}^h,2\nu {\nabla }^s{\mathbf{u}}^h·\mathbf{n}\right)}_{{\Gamma }_{\mathrm{N}}\cap \Gamma }\hfill \\ \hfill & +{\left(q^h,\nabla ·{\mathbf{u}}^h\right)}_{\tilde{\Omega }}+{\mathcal{U}}_{\overline{\Omega }}=0,\forall \{{\mathbf{v}}^h,q^h\}\in {\mathcal{V}}^h\times {\mathcal{Q}}^h.\hfill \end{array}$$

(8)

The operation ${\left(·,·\right)}_{\mathcal{X}}$ refers to the $L^2-$inner product on $\mathcal{X}$= $\tilde{\Omega },{\tilde{\Omega }}_e,\Gamma$, etc., which for two arbitrary vectors $\mathbf{v}$ and $\mathbf{w}$ on the domain $\mathcal{X}$ is given as, ${\left(\mathbf{v},\mathbf{w}\right)}_{\mathcal{X}}={\int }_{\mathcal{X}}\mathbf{v}·\mathbf{w}\mathrm{d}\mathcal{X}$. The last two terms of Equation (8) are, respectively, the continuity constraint and a stabilization term. The stabilization term considered in this work is of the following form:

$${\mathcal{U}}_{\overline{\Omega }}=\sum _{b=1}^{n_{\mathrm{belems}}}{\mathcal{N}}_{{\Gamma }_{\mathrm{b}}\cap {\Gamma }_{\mathrm{D}}}^{\mathrm{Weak}-\mathrm{D}}+{\mathcal{D}}_{{\Gamma }_{\mathrm{out}}\cap \Gamma }^{\mathrm{DDN}}.$$

(9)

This term is the sum of the stabilization of the weakly imposed Dirichlet boundary conditions (see Section 2.2.1), and the directional do-nothing (DDN) condition (see Section 2.2.2).

To discretize the problem, use will be made of both regular B-splines and Truncated Hierarchical B-splines (THB-splines) [45]. THB-splines are particularly useful when performing local mesh refinements, which are crucial in the AFPM application. THB-splines are our local refinement technology of choice because they are available in our software library. We note, however, that alternative refinement technologies such as, PHT-splines [46,47], T-splines [37,48,49], and LR B-splines [50,51], are equally suitable. The B-spline space that is spanned by the above-mentioned basis functions $N_i\left(\xi \right)$ in 1D (with coordinate $\xi$) is denoted by

$$S_{\alpha }^k=\mathrm{span}{\left\{N_i\left(\xi \right)\right\}}_{i=1}^n,$$

(10)

in which k is the spline order and $\alpha$ the regularity. The regularity is related to the degree k and multiplicity $r\le k+1$ by $\alpha =k-r$. In the absence of knot multiplicities within the domain, maximum regularity is achieved, i.e., $\alpha =k-1$, while the boundaries are open with, $\alpha =-1$. A useful property of splines is that their derivatives are also splines, i.e., let $S_{\alpha }^k$ and $S_{\alpha -1}^{k-1}$ be univariate spline spaces, then

$$\left\{\begin{array}{c}\frac{\mathrm{d}}{\mathrm{d}\xi }vs.:vs.\in S_{\alpha }^k\end{array}\right\}=S_{\alpha -1}^{k-1}.$$

(11)

This property is important in ensuring that the velocity space is divergence free [33]. Given that the domain is discretized by rectangular elements, one should specify the element type, i.e., what combination of basis functions are used for the different fields in the different directions. When performing fluid simulations in the context of IGA, Taylor Hood, Nédélec elements of the second family and Raviart–Thomas elements are used [33]. These elements all satisfy the Babuška-Brezzi condition and are therefore stable for the considered mixed formulation. The Raviart–Thomas elements satisfy an additional divergence-free condition, which is ideal for incompressible fluids. The Raviart–Thomas multi-dimensional basis functions are constructed by the tensor-product of multiple univariate basis functions, which for the axisymmetric case yields

$$\begin{array}{cc}\hfill & \mathbf{u}=\left(u_r,u_{\theta },u_z\right) \in {\mathcal{V}}^h:=S_{\alpha ,\alpha -1}^{k,k-1}\times S_{\alpha -1,\alpha -1}^{k-1,k-1}\times S_{\alpha -1,\alpha }^{k-1,k},\hfill \end{array}$$

(12a)

$$\begin{array}{cc}\hfill & p\in {\mathcal{Q}}^h:=S_{\alpha -1,\alpha -1}^{k-1,k-1},\hfill \end{array}$$

(12b)

for, respectively, the discretized velocity and pressure spaces. The elaboration of these spaces is discussed in Remark 1, while various aspects of these spaces are discussed in [33].

Lemma 1.

Given the axisymmetry assumption, the Raviart–Thomas basis of Equation (12) satisfies the divergence-free condition.

Proof .

The general three-dimensional Raviart–Thomas space is equal to

$$\begin{array}{cc}\hfill & {\mathcal{V}}^h:=S_{\alpha ,\alpha -1,\alpha -1}^{k,k-1,k-1}\times S_{\alpha -1,\alpha ,\alpha -1}^{k-1,k,k-1}\times S_{\alpha -1,\alpha -1,\alpha }^{k-1,k-1,k},\hfill \end{array}$$

(13a)

$$\begin{array}{cc}\hfill & {\mathcal{Q}}^h:=S_{\alpha -1,\alpha -1,\alpha -1}^{k-1,k-1,k-1},\hfill \end{array}$$

(13b)

which ensures the Babuška–Brezzi condition. The axisymmetry definition equals

$$\frac{\partial ·}{\partial \theta }=0,$$

(14)

which states that the solution is independent of the position in the $\theta$-direction (see Figure 4). Inserting this axisymmetry condition in Equation (13) yields the Raviart–Thomas basis of Equation (12). □

2.2.1. Weakly Imposed Dirichlet Boundary Condition

The Dirichlet boundary conditions (BC) are weakly imposed instead of enforcing them through constraints on the function space.The most prominent advantage of this weak imposition is that it avoids oscillatory behavior near boundaries when higher-order basis functions are used [52]. The weak imposition of boundary conditions was first derived by Nitsche for the Poisson Equation [53]. Application to the Navier–Stokes equation is discussed by several authors [52,54], of which the terms and constants proposed by Hughes and Bazilevs [29] are considered in the remainder of this manuscript.

The idea of Nitsche’s method is to add additional terms to the variational equations to enforce Dirichlet BCs as Euler–Lagrange conditions. The Nitsche terms to be added to the variational form of the Navier–Stokes Equation (3) are defined as

$$\begin{array}{cc}\hfill {\mathcal{N}}_{{\Gamma }_{\mathrm{b}}\cap {\Gamma }_{\mathrm{D}}}^{\mathrm{Weak}-\mathrm{D}}=& -{\left({\mathbf{u}}^h,2\nu {\nabla }^s{\mathbf{v}}^h·\mathbf{n}\right)}_{{\Gamma }_b\cap {\Gamma }_{\mathrm{D}}}\hfill \\ & -{\left(\gamma 2\nu {\nabla }^s{\mathbf{v}}^h·\mathbf{n},{\mathbf{u}}^h-{\mathbf{u}}_{\mathrm{D}}^h\right)}_{{\Gamma }_b\cap {\Gamma }_{\mathrm{D}}}\hfill \\ & +{\left(\frac{C_b^I\nu }{h_b}{\mathbf{v}}^h,{\mathbf{u}}^h-{\mathbf{u}}_{\mathrm{D}}^h\right)}_{{\Gamma }_b\cap {\Gamma }_{\mathrm{D}}},\hfill \end{array}$$

(15)

in which $h_b$ is the element normal length at the boundary, ${\Gamma }_{\mathrm{D}}$ the boundary at which the velocity is assigned a Dirichlet condition and ${\mathbf{u}}_{\mathrm{D}}^h$ is the desired velocity vector at the wall. Note that Equation (15) is only defined at the element level and should still be summed over all boundary elements in accordance with Equation (9). Following the arguments in Ref. [29], the constant $\gamma$ is taken equal to 1. The constant $C_b^I$ should be sufficiently high and dependent on the order of interpolation and element type used. In the upcoming analyses, this value is taken equal to $C_b^I=5p^{-2}$, which is based on Ref. [55] and further adjusted empirically as it showed satisfactory results for the considered test cases. Note that Nitsche’s method is consistent with the original problem, since the latter two terms in Equation (15) vanish for ${\mathbf{u}}^h={\mathbf{u}}_{\mathrm{D}}^h$ at ${\Gamma }_{\mathrm{D}}$. An in-depth derivation of Nitsche’s method in combination with an analysis regarding convergence and stability can be found in Refs. [52,54,56].

2.2.2. Directional Do-Nothing Boundary Condition

In many flow problems, including the one considered here, an infinite ambient domain is bounded for computational reasons. In the case that such a boundary has neither a Dirichlet nor a non-homogeneous Neumann BC assigned to it, it is considered as an open or do-nothing boundary. Frequently, such a boundary corresponds to a pure outflow and does not require any further treatment. However, in the case of a swirling or recirculating flow, a pure outflow cannot be guaranteed as possible inflow might occur. Special treatment is then required as such inflow may lead to numerical stability issues. To circumvent such issues, we use the modification of the do-nothing condition by Braack and Mucha [57] defined as

$$\begin{array}{cc}\hfill {\mathcal{D}}_{{\Gamma }_{\mathrm{open}}\cap \Gamma }^{\mathrm{DDN}}& =-{\left({\mathbf{v}}^h,\left[2\nu {\nabla }^s{\mathbf{u}}^h-p^h\mathbf{I}\right]·\mathbf{n}\right)}_{{\Gamma }_{\mathrm{out}}}\hfill \\ & =-\frac{1}{2}{\int }_{{\Gamma }_{\mathrm{out}}}{\left({\mathbf{u}}^h·\mathbf{n}\right)}_{\mathrm{DDN}}{\mathbf{u}}^h·{\mathbf{v}}^h\mathrm{d}\Gamma ,\hfill \end{array}$$

(16)

where

$${\left({\mathbf{u}}^h·\mathbf{n}\right)}_{\mathrm{DDN}}:=\left\{\begin{array}{c}0\mathrm{for}{\mathbf{u}}^h·\mathbf{n}\ge 0,\\ {\mathbf{u}}^h·\mathbf{n}\mathrm{for}{\mathbf{u}}^h·\mathbf{n}<0.\end{array}\right.$$

(17)

Equation (16) satisfies the classical do-nothing condition when outflow occurs, but vanishes when inflow is present. Due to it being direction dependent, this condition is referred to as the directional-do-nothing (DDN) BC. The dependence of this term on the flow direction ensures stability [57].

2.3. Validation

In this section, the weak formulation presented above is validated using three different flow problems for which analytical solutions are available. To be representative for the AFPM machine application, all problems pertain to an axisymmetrically represented rotating disk. It is also chosen to strongly impose the normal velocity component on the wall boundary conditions in accordance with Ref. [29].

2.3.1. Enclosed Rotating Disk

The first flow problem is comprised of a rotor-stator configuration with a shroud, i.e., a closed-end configuration (see Figure 5). The domain consists of a symmetry boundary, which only requires a no-penetration condition, while the remaining wall boundaries require additional constraints. The top wall has a prescribed rotational speed, which is selected such that the Reynolds number based on the rotor radius is $\mathrm{Re}<10$. The model-specific parameters are given in

Table 1. Parameters used for the enclosed rotating disk test-case.

$\mathbf{\Omega }$ [rad/s]	$\mathit{\nu }$[m${}^2$/s]	Domain Size $\mathit{R}\times \mathit{H}$ [m]	Element nr. $\mathit{R}\times \mathit{H}$ [-]	Spline Order [-]
1.0	1.0	$1.0\times 0.5$	$30\times 15$	2

. The low Reynolds number makes it possible to neglect the convection. One can then derive an analytical 2D axisymmetric steady-state solution [58], which is provided in Appendix A.

Due to the absence of out- and inflow boundaries and the low Reynolds number, the Directional Do-Nothing term has been neglected in the numerical model. Comparison between the analytical solution and the numerical model in Figure 5b shows excellent agreement. This problem demonstrates the effectiveness of the weakly imposed Dirichlet boundary conditions. The effect of the weak imposition is particularly visible in the top right corner, where the rotor and wall boundary meet. At this point, the velocity is ill-defined as it cannot satisfy both the rotor and wall boundary condition. Strongly imposed Dirichlet conditions, in such cases, cause oscillations, while the weak variant maintains stability.

2.3.2. Von Kármán Swirling Flow

The swirling flow problem is a well-known problem in the context of rotating disk flows. This problem was first proposed and (analytically) solved by von Kármán [59]. Von Kármán considered the flow induced by an infinite rotating disk where the far-field fluid is at rest. By means of the self-similarity principle, he was able to reduce the Navier–Stokes equations to a set of nonlinear ordinary differential equations. The derived solution can be expressed in four dimensionless quantities

$$u_r=r\Omega F\left(\zeta \right),u_{\theta }=r\Omega G\left(\zeta \right),u_z=\sqrt{r\Omega }H\left(\zeta \right),$$

(18)

with

$$\zeta =z\sqrt{\Omega /\nu }.$$

(19)

These quantities, respectively, represent the radial, tangential and axial velocity as a function of the dimensionless axial height $\zeta$. $\Omega$ is the rotational speed of the disk, $\nu$ the kinematic viscosity and r and z are the radial and axial coordinates, respectively. The domain used for our numerical analysis is illustrated in Figure 6a, while the model-specific parameters are listed in

Table 2. Parameters used for the von Kármán swirling flow test-case.

Time-Step [s]	$\mathbf{\Omega }$ [rad/s]	$\mathit{\nu }$ [m${}^2$/s]	Domain Size $\mathit{R}\times \mathit{H}$ [m]	Element nr. $\mathit{R}\times \mathit{H}$ [-]	Spline Order [-]
0.02	1.0	0.2	$12.0\times 12.0$	$10\times 30$	2

.

It is important to note that the von Kármán solution is only valid in an infinite domain, which is physically impossible to simulate. In order to minimize the boundary influence in our finite-domain simulations, the problem is solved transiently. The model is considered convergent when the boundary layers are fully developed, i.e., they do not change substantially with increasing time. For the considered parameters, the Reynolds number is equal to 720 at the outer radius of the disk. The numerical solution is monitored at a fixed radial position of $r/r_{\max }=0.4$, where the local Reynolds number equals 125. Note that the right boundary is open (outflow), without constraints. Hence, the directional do-nothing condition is used in this simulation. The numerical results shown in Figure 6b are in good agreement with the analytical solution [59].

2.3.3. Infinite Rotor-Stator Configuration

The final validation case is effectively a combination of the former two. It consists of a rotating disk with a rotational speed of $\Omega =1$ and a stationary disk, separated by a disk spacing d as illustrated in Figure 7a. The radii of both the rotor and the stator are assumed to be infinite, meaning that no complex outflow behavior is present. The reference solution is derived in Ref. [60] using a power series expansion. The dimensionless numbers required for the solution are

$$\widehat{z}=\frac{z}{d},\widehat{\mathbf{u}}=\frac{\mathbf{u}}{\Omega d}=\left[r{H}^{\left(1\right)}\left(\widehat{z}\right),rG\left(\widehat{z}\right),-2H\left(\widehat{z}\right)\right],$$

(20)

and

$${\mathrm{Re}}_d=\frac{\Omega d^2}{\nu }.$$

(21)

The results are generated using the parameters listed in

Table 3. Parameters used for the infinite rotor-stator configuration test-case.

$\mathbf{\Omega }$ [rad/s]	$\mathit{\nu }$[m${}^2$/s]	Domain Size $\mathit{R}\times \mathit{H}$ [m]	Element nr. $\mathit{R}\times \mathit{H}$ [-]	Spline Order [-]
1.0	1.0	$2.0\times 1.0$	$15\times 15$	2

, and are presented in Figure 7d, which correspond to ${\mathrm{Re}}_d=5$.

The radial velocity in Figure 7b clearly shows the von Kármán boundary layer at the rotor and the Bödewadt boundary layer at the stator. The outward moving fluid at the rotor is accompanied with an inflow at the stator, satisfying the mass conservation. Due to the low Reynolds number, the tangential velocity in Figure 7c shows an expected linear torsional Couette flow, in agreement with the literature [10]. The remaining axial component shown in Figure 7d is small but non-zero. Increasing the Reynolds number eventually results in a non-viscous core, as explained in the introduction and shown in Ref. [60]. Reproducing this non-viscous core is, however, difficult to achieve with the current model, as stable simulations at these high Reynolds number flows are not obtained with the considered formulation. The laminar results show a very good agreement with the power series solution, however, validating the model for the low Reynolds number regime considered in this work.

3. Thermal Model

Isogeometric analysis is also used to solve the thermal problem on an axisymmetric domain. In our solution procedure, the thermal model is decoupled from the flow model, i.e., the flow field is isothermally computed after which it serves as an input for the temperature field calculation. This decoupling assumption is valid as long as the temperature differences in the domain do not alter the density and viscosity significantly (see Appendix B).

3.1. Strong Form

The thermal energy equation is a simplification of the more general (total) energy balance, which also takes potential and kinetic energy into account. While the potential and kinetic energy are very important when dealing with compressible flows at high speeds (Mach number, $\mathrm{Ma}>0.3$), they can be neglected in the current application. The reason for this is the order of magnitude difference between the thermal energy, $c_p\Delta T\sim \mathcal{O}\left({10}^4\right)$ [J/kg], and the kinetic energy, $0.5v^2\sim \mathcal{O}\left(1\right)-\mathcal{O}\left(10\right)$ [J/kg], where $c_p$ is the specific heat capacity. Hence, kinetic energy effects can safely be disregarded. The energy balance then reduces to a convection-diffusion equation, which in a steady-state takes the following form

$$\nabla ·\left(\rho c_p\mathbf{u}T\right)=\nabla ·\left(k\nabla T\right)+Q_S,$$

(22)

with the boundary conditions

$$T=g\mathrm{at}{\Gamma }_{\mathrm{D}}\mathrm{and}\mathbf{n}·\nabla T=0\mathrm{at}{\Gamma }_{\mathrm{N}}.$$

(23)

In these expressions, T [K] is the unknown temperature and $\mathbf{u}$ [m/s] the velocity vector, which is pre-computed, as discussed in Section 2. The boundary conditions follow the definitions in Figure 4. The volumetric heat source or sink term $Q_S$ [W/m${}^3$] encompasses the heat generated by the eddy current losses in the core. The thermal conductivity k [W/mK] is a material (fluid or solid) property depending on the position inside the domain. In line with the assumption that temperature differences are relatively small, it is assumed that the conductivity of both the fluid and the solid material is independent of the temperature.

3.2. Variational Form

As discussed in the previous section and illustrated in Appendix B, the velocity, density, and viscosity coupling is ignored on account of the assumed small variation in the temperature. Consequently, the system of equations is unidirectionally coupled through the velocity field. This coupled system is here solved using a segregated approach, which solves the equations one after another. Although this might require multiple iterations before the system of equations converges, it requires far less memory and could, therefore, be faster compared to the fully coupled approach.

To derive the variational formulation for the convection-diffusion problem in Equations (22) and (23), we define the discretized trial and test space as

$${\mathcal{T}}^h:=\{T^h\in H^1\left(\tilde{\Omega }\right)|T^h=T_{\mathrm{D}}\mathrm{on}{\Gamma }_D\},$$

(24)

$${\mathcal{W}}^h:=\{w^h\in H^1\left(\tilde{\Omega }\right)|w^h=0\mathrm{on}{\Gamma }_D\},$$

(25)

in which the subscript D refers to the Dirichlet boundaries. For an incompressible flow, the discretized variational formulation on an arbitrary domain $\tilde{\Omega }$ then reads:

Find $T^h\in {\mathcal{T}}^h$ such that:

$$\begin{array}{cc}\hfill & {\left(w^h,\left[c_p\rho \mathbf{u}\right]·\nabla T^h\right)}_{\tilde{\Omega }}+{\left(\nabla w^h,k\nabla T^h\right)}_{\tilde{\Omega }}\hfill \\ & -{\left(w^h,Q_S^h\right)}_{\tilde{\Omega }}=0,\forall w^h\in {\mathcal{W}}^h.\hfill \end{array}$$

(26)

For the AFPM application, there is no heat flux through the domain boundaries. Hence, merely Dirichlet and homogeneous Neumann boundary conditions are applied. Stabilization techniques such as the Streamline-Upwind Petrov-Galerkin (SUPG) method are not required in this work as the considered Peclet number is relatively low. As for the flow model discussed in Section 2, the thermal problem is discretized using regular and truncated hierarchical B-splines.

In this work we consider a conjugate heat transfer problem, i.e., the temperature distribution in the solid disks is modeled, using Equation (22), while ignoring convection. If radition becomes a dominant heat transfer mechanism, this would require consideration of a Surface-to-Surface, P-1, or Rosseland model, which simulates the various radiative effects between the different parts in a domain. Extension of our modeling approach with such radiative effects is considered a topic of future study.

3.3. Nusselt Number

As discussed briefly in Section 1.2, the dimensionless Nusselt number is often used to indicate the relation between convective and conductive heat transfer in a system. The Nusselt number is especially useful when comparing the effectiveness of geometrical improvements, such as cooling vanes, as it is solely evaluated at the surface of an object. For a disk-shaped geometry, the Nusselt number is defined as in Equation (2). Given that the convective heat transfer coefficient is dependent on the velocity, it can be related to the Reynolds number. Although the relation between the velocity and the heat transfer coefficient is problem dependent, an increase in Reynolds number will also increase the Nusselt number, due to the increase in convective heat transfer. Since the definition in Equation (2) is not practical when dealing with complex geometries, an alternative averaged Nusselt number can be defined as

$$\overline{\mathrm{Nu}}=\frac{q_{\mathrm{av}}R}{k{\left(T_{\mathrm{wall}}-T_0\right)}_{\mathrm{av}}}=\frac{\int r\nabla T·\mathbf{n}\mathrm{d}r}{\int r\left(T_{\mathrm{wall}}-T_0\right)\mathrm{d}r}R,$$

(27)

in which $q_{\mathrm{av}}$ [W/m${}^2$] is the average heat flux through the heated surface and $T_{\mathrm{wall}}$ [K] the temperature of that considered surface. The Nusselt number (27) in essence divides the average heat transfer in the fluid by its average wall temperature. Following Howey [5], the reference temperature $T_0$ should be considered as a function of the bulk flow temperature. Based on the energy balance across the rotor-stator gap, one obtains

$$T_0\left(r\right)=\frac{1}{\dot{m}}\int \int \rho u\left(r,z\right)T\left(r,z\right)\mathrm{d}\theta \mathrm{d}z,$$

(28)

with $\dot{m}$ [kg/s] the mass flow, u [m/s] the velocity parallel to the two disks and T [K] the temperature. This definition is experimentally highly inconvenient for an AFPM application, since one requires the flux across an axial gap at which recirculation might occur [5]. To avoid these complications, in the remainder of this work the reference temperature is taken equal to the ambient temperature, although this might result in the Nusselt number to become dependent on the fluid temperature.

3.4. Validation

The validation of both the temperature field and the Nusselt number is based on the von Kármán swirling flow discussed in Section 2.3.2. This validation problem consists of a rotating disk for which the flow far away is at rest. Application of the thermal energy Equation (26) enables the calculation of the temperature distribution and consequently the average Nusselt number (27). The domain is presented in Figure 8a, while the model-specific parameters are given in

Table 4. Parameters used for the von Kármán rotating disk heat transfer test-case.

$\mathbf{\Omega }$ [rad/s]	$\mathit{\nu }$ [m${}^2$/s]	Domain Size $\mathit{R}\times \mathit{H}$ [m]	Element nr. $\mathit{R}\times \mathit{H}$ [-]	Spline Order [-]	Prandtl Number [-]	$\mathit{\beta }$ [-]	${\mathit{T}}_0$ [${}^{\circ }$C]
1.0	0.2	$12.0\times 12.0$	$10\times 30$	2	0.72	100.0	20.0

. The result of the temperature distribution is shown in Figure 8b, where use is made of the following dimensionless parameters:

$$\Theta =\frac{T-T_0}{T_{\mathrm{disk}}-T_0},\zeta =z\sqrt{\Omega /\nu },\mathrm{Pr}=\frac{\mu c_p}{k}.$$

(29)

The Prandtl number, Pr, is a fluid-specific value that relates the momentum or viscous diffusivity to the thermal diffusivity. The numerical results of Figure 8b, show excellent agreement with the analytical power-law solution by Shevchuk [61] and the experimental results of Elkins [62].

The Nusselt number validation is based on the relation provided by [63], which uses a power-law relation such that

$$\overline{\mathrm{Nu}}=a{\mathrm{Re}}^{0.5},$$

(30)

in which a is a constant dependent on the Prandtl number, temperature distribution, and flow regime. Exact values of a are provided in Ref. [61], based on the similarity solution. Suppose a radial disk temperature distribution is prescribed using the power-law

$$T_{\mathrm{disk}}\left(r\right)=T_0+\beta r^n,$$

(31)

where $T_0$ [K] is the reference temperature and $\beta$ [-] an arbitrary non-zero constant (see Table 4). For different values of n, one can then solve the temperature distribution and compute the Nusselt number. The result of the model in comparison with the exact similarity solution is provided in Figure 9. This figure conveys that the obtained simulation results are in very good agreement with the reference result, and demonstrates the strong dependence of the Nusselt number on the radial temperature distribution. The slight discrepancy at higher values of n is related to the computed velocity field, which is influenced by the outflow boundary. An increase in n shifts the disk temperature distribution towards the periphery, where the velocity field is affected by its boundary (i.e., circulation occurs). This ultimately influences the average Nusselt number, as shown in Figure 9.

4. AFPM Simulation Setup

The simulations considered in the remainder of this work are based on the experimental setup, illustrated in Figure 10, which differs slightly from a typical AFPM machine.

The current-carrying coils, located at the stator, are removed for experimental purposes. An external motor is then required to drive the rotor, illustrated in Figure 10b. As the magnet array rotates, eddy currents are induced in the iron core, which generate a damping torque and produce heat. By varying the rotation speed and gap-length, the influence of the core geometry can be investigated using the experimental setup. The flow and temperature field are, however, not influenced by the change in electromagnetic fields compared to the actual AFPM machine. Note that the illustrated rotor and stator in Figure 10a are both complex 3D periodic geometries, i.e., a specific pie-shaped part of the disk geometry repeats itself in the $\theta -$direction. This periodicity is not accounted for in our simulations as a consequence of the axisymmetry assumption, hence the angular slits in the core and between the magnets as seen in Figure 10b are not modeled geometrically. Although it would be favorable to include such three-dimensional features, this is impractical from a computational effort point of view. Based on results reported in the literature, it is expected that the axisymmetric model will result in similar flow and heat transfer structures compared to the three-dimensional case [4,64].

The geometry is furthermore simplified inside the rotation plane itself. The magnet array, back-iron, and core are all modeled inside the rotor and stator. This change ultimately influences the flow behavior but is considered practical from a mesh generation point of view. The experimental setup also has optional circular air-inlets, periodically positioned at the stator disk with a specific diameter. When assuming axisymmetry, these holes will result in an extruded ring. This problem simplification should be considered when analyzing the results, as a larger inflow area is created. Both disks are assumed to have an equal radius of $R=0.17$ m. The effect of buoyancy, i.e., natural convection, is neglected because of two reasons. First, the temperature difference in combination with the occurring convection is expected to be sufficient to mitigate most effects of buoyancy. Second, the axis of rotation of the experimental setup is positioned perpendicular to the gravitational direction. This limits the model to capture any gravitational buoyancy as a result of the axisymmetry assumption.

4.1. Numerical Domain and Mesh

The numerical domain is based on an open-end rotor-stator configuration, in combination with adequate boundary conditions (BC). An illustration is provided in Figure 11, where different materials correspond to different colors. Besides the rotor and the stator, a substantial part of the ambient environment is modeled. This is required for open-end rotor-stator configurations [17,24].

A distinction is made between the flow and temperature domain. Both are equal in size, but the temperature mesh includes the solid material, while the flow domain excludes the rotor and stator parts, i.e., the Navier–Stokes and continuity equations are not required to be solved here. The boundaries of the rotor and the stator are then subject to a prescribed wall velocity, ${\mathbf{u}}_{rotor}$ and ${\mathbf{u}}_{stator}$, respectively. The stator core experiences eddy losses in the form of heat. Hence, an external heat source $Q_{core}$ is assigned to the core domain, which is indicated in orange in Figure 11. Because the current multi-physics model is decoupled from the electro-magnetic model used in Ref. [2], the heat source is chosen to be uniformly distributed, which simplifies the upcoming analyses. In reality, however, this heat source should be non-uniform due to eddy current losses. More specifically, a maximum is expected in the center of the core, influencing the temperature distribution. The magnitude of the uniform heat source is chosen such that the average temperature of the surface at the gap side of the core equals a prescribed temperature

$${\overline{T}}_{\mathrm{c}}=\frac{1}{A_{\mathrm{c}}}\underset{0}{\overset{2\pi }{\int }}\underset{R_{i,c}}{\overset{R_{o,c}}{\int }}r{T}_{\mathrm{c}}\left(r\right)\mathrm{d}r\mathrm{d}\theta ,$$

(32)

with

$$A_{\mathrm{c}}=\pi \left(R_{o,c}^2-R_{i,core}^2\right),$$

(33)

the area of the core-surface in the tangential and radial plane. Note that ${\overline{T}}_{\mathrm{c}}$ is solely defined at the surface that is in contact with the gap. The remaining thermal BCs consist of a fixed ambient temperature, $T_0$, at the left and right boundaries of the domain and an open BC at the top boundary. These boundary conditions are set in accordance with the observed fluid flow, which acts as an inlet at the left and right boundaries while being an outlet at the top. The ambient temperature and average core surface temperature are set to $20{}^{\circ }$C and $120{}^{\circ }$C, respectively. The axis of rotation uses a no-penetration as well as an open BC, similar to the validation test-cases discussed above. The fluid properties are based on air at $80{}^{\circ }$C and are assumed constant. The properties used in the simulations are listed in

Table 5. Parameters used for the AFPM simulations.

Disk Radius [m]	$\mathbf{\Omega }$ [rad/s]	$\mathit{\nu }$ [m${}^2$/s]	$\mathit{\rho }$ [kg/m${}^3$]	${\mathit{c}}_{\mathit{p}}$ (air) [J/kgK]	Domain Size $\mathit{R}\times \mathit{H}$ [m]	Spline Degree [-]
0.17	0.0–6.0	$2.097\times {10}^{-5}$	1.0	1008.0	$2.5\times 5$	2.0

. Additional information regarding the thermal conductivities is given in Appendix C.

The domain size is selected such that the ambient boundary conditions do not influence the flow behavior near the disks. The element size is increased toward the far-field boundaries to reduce the total number of elements in the mesh. The possible reduction in accuracy near the boundaries does not significantly affect the region of interest.

Typical results for both the velocity and temperature fields are provided in Figure 12a,b, respectively. Some small velocity oscillations are observed above the rotor (see Figure 12a) which are related to the element size and Reynolds number. Further reduction of the element size would remove such oscillations. Provided that the oscillations remain small and outside the region of interest, they are expected not to significantly influence the quantities of interest.

4.2. Mesh Convergence

The suitability of the mesh is studied by monitoring the quantities of interest under mesh refinement. These quantities correspond to the radial velocity component inside the gap and the average Nusselt number across the heated core (27). Note that the radial velocity component will contribute the most in terms of cooling, provided that axisymmetry is assumed and the axial component is relatively small. For the mesh convergence study presented here, the geometry was chosen to have a gap-length of 5 mm, without the addition of an air-inlet. Boundary conditions correspond to a maximum rotation speed of $4.0$ rad/s, yielding a Reynolds number of $\mathrm{Re}=5513$. Since we are mainly interested in the heat transfer inside the gap, this region is refined. By dyadically refining the gap, the results of Figure 13 and

Table 6. Number of elements and dofs for the mesh convergence study for the fluid domain. The reported number of degrees of freedom pertains to quadratic spline basis functions.

	${\mathcal{M}}_1$	${\mathcal{M}}_2$	${\mathcal{M}}_3$	${\mathcal{M}}_4$
$n_{el}$ [-]	20,256	21,216	29,280	53,640
$n_{dof}$ [-]	82,252	86,484	119,388	217,914

were obtained, with a spline order of 2. Note that since only the gap is refined, the number of elements (and degrees of freedom) does not increase by a factor of 4 as one would expect. In particular for the coarsest meshes, the number of elements inside the gap is small compared to the total number of elements. The substantial increase in number of elements on the finest mesh, ${\mathcal{M}}_4$, is caused by the fact that it required the surrounding environment (not only the gap) to be refined as well, in order to obtain a gradual mesh transition.

A detailed visualization of the ${\mathcal{M}}_3$ fluid mesh is provided in Figure 14, clearly showing the difference in mesh size across the domain as well as the element stretching near the boundaries. Computation times ranged from 6 up to 8 h excluding mesh creation, on a Linux cluster with 1TB available RAM and 4 CPUs. The system matrix assembly is performed in parallel, while the solving procedure is performed in serial. The latter makes use of the Math Kernel Library (MKL) [65], which includes an efficient serial direct solver. Computation times were mainly dominated by the temperature calculation due to the large number of dofs inside the solids that were required to capture the thin core segment. The computation times consisted of roughly 20% pre-processing, 75% calculation, and 5% post-processing. These times will inevitably increase when higher rotation speeds are required, as multiple Newton increments will then be required. The results in Figure 13 show clearly that the third mesh, ${\mathcal{M}}_3$ strikes a good balance between accuracy and computation time. Therefore, this mesh resolution will be used inside the gap for the simulations in the remainder.

5. AFPM Design Analyses

Using our computational model, the influence of the rotation speed and gap-length on the heat transfer inside the gap is studied. Two AFPM configurations are considered, viz., one with and one without an air-inlet. The observed thermal behavior is explained by means of the flow structures present in the AFPM machine, leading to design recommendations pertaining to the cooling of the machine.

The gap-length d shown in Figure 3 is varied within a specific range to analyze the effects on the heat transfer. The range is based on the experimental setup which allows for a gap-length of $1.0$ mm at its minimum and $20.0$ mm at its maximum. Relatively low rotation speeds will be monitored in order to remain within the laminar flow regime to which our model is applicable.

Solving the non-linear Navier–Stokes equation on the considered complex geometry required an incremental-iterative Newton–Raphson method. This method incrementally increases the rotation speed, while using the previously determined field as an initial guess for the next increment. This ensured the solver to reach a global minimum, up until convection-related instabilities would occur. The provided results are limited to a maximum Reynolds number of $\mathrm{Re}=8269$. Simulating higher Reynolds number flows requires extension of the computational model with additional stabilization methods, such as the Streamline-Upwind Petrov-Galerkin (SUPG) method, time-dependent solvers, and turbulence models. The extension to higher Reynolds number flows is of interest for the considered application, but is considered beyond the scope of the work presented here.

The results for the heat transfer across the core surface inside the gap are provided in Figure 15, showing the differences between the configurations with and without air-inlet. Note that the Nusselt number is relatively low for all considered cases in Figure 15. There are two reasons for this: (I) The Reynolds number is substantially lower than some reported in the literature (which corresponds to $\mathrm{Re}\sim \mathcal{O}\left({10}^5\right)$ resulting in $\overline{\mathrm{Nu}}\sim \mathcal{O}\left({10}^2\right)-\mathcal{O}\left({10}^3\right)$). Such high Reynolds numbers can only be obtained when using a suitable turbulence model, which is not discussed here. (II) The Nusselt number defined in Equation (27) uses an incorrect reference temperature definition by taking it equal to the ambient temperature. This is correct for geometrically simple problems, but becomes problematic in the setting considered here. The definition in Equation (28) would correct this, but this definition is deemed impractical should experiments be conducted as it requires the temperature across the gap to be known. By using the ambient temperature, an over-estimated temperature difference of $\Delta T=100{}^{\circ }$C is found, while the correct definition equals roughly $\Delta T=$ 30–40 ${}^{\circ }$C. This is because the average bulk temperature inside the gap is larger than the ambient reference temperature. Such an increase in the temperature difference ultimately decreases the Nusselt number in accordance with Equation (27).

Figure 15 shows a steep increase in heat transfer when the gap-length ratio, G (see Equation (1)), decreases. Larger values of G on the contrary, do not show such steep variations. The non-uniform temperature boundary condition employed here leads to different behavior compared to Yuan et al. [24], who observe a decrease in Nusselt number upon decreasing the aspect ratio at low gap ratios. More specifically, they observe an optimum (maximum) Nusselt number at a specific aspect ratio. Because of the uniform temperature across the entire stator as considered by Yuan et al., the entering fluid is continuously heated, without the ability to dissipate any thermal energy. This leads to supposedly low Nusselt numbers at small gap sizes. To clarify this further, both the flow and the temperature profiles are analyzed in more detail below. Use will be made of a fixed Reynolds number ($\mathrm{Re}=5513$, corresponding to the blue curve in Figure 15), while the gap-lengths are varied.

5.1. Flow Field

Figure 16 and Figure 17 show the flow field present inside the AFPM machine. Use is made of the Line Integral Convolution (LIC) visualization technique available in Paraview [66]. The visualized gap is a sub-domain of the entire numerical AFPM domain, as shown in Figure 12a by the box, $\mathbf{B}$. Note that the flow direction is dictated by the fact that the flow at the rotor side is always outwards, i.e., it flows from the left to the right. When observing the configuration without an air-inlet in Figure 16, a similar flow trend is seen for all aspect ratios, G. At the rotor side, the air is pumped outwards while an equal amount is sucked inwards at the stator side. As a result, a circulation zone appears near the periphery, as also observed by Yuan et al. and Debuchy et al. [17,24]. Both research results also showed that the center of the circulation zone is slightly more positioned towards the rotor disk, which is also seen in Figure 17 for larger aspect ratios.

The configuration with air-inlet shows more complex flow behavior (see Figure 17). The smallest gap-length ratio shows a dominated outflow inside the gap, caused by the air-inlet, which results in a positive (outward) flow direction at the stator side. However, a small amount of fluid is still sucked inward at the periphery, creating a minor circulation zone as also observed in Figure 16. In contrast to Figure 16, the circulation is mainly located at the stator, rather than at the rotor disk. The size of this circulation zone increases with increasing gap-length aspect ratio, causing the flow direction to change at the stator. If the gap-length ratio is sufficiently large, the circulation zone effectively approaches the flow field as also seen for the configuration without air-inlet. This circulation zone has also been observed by Howey [5] and Soo [67], who investigated similar configurations. It is also present inside a closed-end rotor-stator configuration with superposed throughflow as is shown by Poncet [15]. Superposed here refers to a forced inflow of fluid at the air-inlet, in contrast to the open throughflow present inside the current AFPM machine. This circulation zone is highly influenced by the inflow velocity at the air-inlet, but also by the gap-length and the Reynolds number. Specific combinations of these three parameters will result in different flow fields and heat transfer characteristics. The effect of varying the gap-length ratio is already presented in Figure 17. The influence of the Reynolds number for a specific gap-length has also been investigated. The results showed that an increase in Reynolds number causes the circulation zone to shift towards the axis of rotation. However, this effect is minor, and hence the results are not presented in detail here. Note that Figure 17 also shows a circulation zone near the axis of rotation upon decreasing the aspect ratio. This effect does, however, not influence the flow near the core.

The velocity profiles in tangential and radial direction are visualized in Figure 18 for three different radial positions $r^{*}=r/R$, viz., (a) $0.52$, (b) $0.73$, and (c) $0.93$. These correspond to the inner, middle, and outer location of the core in the radial direction, as also visualized in Figure 16 and Figure 17. Since the axial velocity is relatively low in magnitude and therefore has limited influence on heat transfer, it is not shown in these figures. A clear distinction can be made between the various observed flow profiles. As explained in the introduction, two kinds of flow types can be expected, namely Batchelor flow and Stewartson flow. The difference between both is that the former has a boundary layer on both the rotor and the stator in the tangential direction, while Stewartson has a boundary layer at the rotor only. Figure 18 clearly shows no boundary layer at the stator disk for the normalized tangential velocity $u_{\theta }^{*}$. Hence, the flow present in both configurations is of the Stewartson type, for all gap-length ratios. It is also observed that the smallest gap-length ratio exhibits a torsional Couette flow. It is expected that a Batchelor flow emerges when a sufficiently high Reynolds number is achieved with a relatively high gap-length, but this regime is beyond the reach of the current model.

The influence of an air-inlet is mainly visible in the radial velocity component, $u_r^{*}$. More specifically, the radial velocity for the gap-length ratio $\mathrm{G}=0.03$, i.e., the blue dotted line in Figure 18a, corresponds to a Hagen–Poiseuille velocity profile. The influence of the air-inlet starts to disappear when approaching the periphery at $r^{*}=1$ due to the above-mentioned circulation zone.

5.2. Temperature Field

The temperature distribution along the stator wall and inside the gap is illustrated in Figure 19 and Figure 20. The radial distribution shown in Figure 19 provides insight into the cooling of the stator, while Figure 20 provides a comparison of the temperature gradients across the gap.

A drop in temperature is noticed when approaching the axis of rotation, $r^{*}\to 0$. This is the result of heat transfer taking place at surfaces other than the core surface. More specifically, heat is dissipated to the ambient domain at the remaining stator sides. It is therefore important to understand that the core temperature distribution, which is chosen on average to be equal to $\overline{T}=120$ ${}^{\circ }$C, is not only cooled by the flow at the core surface. For example, when considering the air-inlet configuration, a major improvement in the overall stator cooling is observed when $r^{*}<0.5$, compared to the counterpart without air-inlet. At first sight, this might contradict the heat transfer results of Figure 15, which shows a decrease in heat transfer when using an air-inlet in the majority of cases. It is, however, important to realize that the Nusselt number in Figure 15 quantifies the heat transfer due to the flow present at the core-surface only. The heat transfer at the core surface is only a part of the overall cooling of the stator disk. This is observed in Figure A2, Appendix D, where a much higher core heat source was required to obtain the same average core surface temperature. Even though equal-sized air-inlets are used for all considered cases, a difference in cooling is observed, with $\mathrm{G}=0.06$ being the most efficient as it reaches the lowest temperature across the entire surface. These differences are attributed to a difference in velocity magnitude present for each configuration. More specifically, there is an optimum aspect ratio, where the inflow velocity at the air-inlet is at its maximum for a specific Reynolds number, which in this case corresponds to $\mathrm{G}=0.06$.

Regarding the heat transfer at the core surface of the configuration without air-inlet, it is noticed that both $\mathrm{G}=0.06$ and $0.12$ yield similar results, as they are mainly cooled at the periphery due to the ingress of cooler ambient air. The smaller aspect ratio, $\mathrm{G}=0.03$, in contrast, shows a different behavior, with the maximum temperature in the middle of the core. This indicates that cooling also occurs at the inner radius of the core, which is a consequence of the relatively thin boundary layer and the increase in surface-area-to-volume (SA-V) ratio due to the decrease of G. Both aspects are beneficial in terms of cooling, as the heat transfer rate increases with increasing SA-V ratio. A lower cooling effect is seen at the periphery, as the complete boundary layer heats up faster compared to the larger gap-length ratios. Toward the inner side of the core, the temperature decreases again as it interacts with the flow at the rotor, which has been cooled when passing the inner part of the rotor and stator.

The configuration with air-inlet behaves differently, as the flow becomes mainly outward at small gap-length ratios. The absence of a circulation zone in combination with the increased SA-V ratio, results in an increase in heat transfer, compared to larger gap-length ratios. This explains why an increase in heat transfer was observed in Figure 9 upon decreasing the gap-length ratio. The difference with the research of Yuan et al. [24], who observed a decrease in heat transfer, is related to the stator boundary conditions. As discussed above, the uniform temperature condition in Ref. [24] causes the fluid to continuously rise in temperature upon entering the system. This would inevitably result in a decrease in heat transfer, which explains the difference with our simulations. It should be noted that both results are in entirely different flow regimes. Extrapolation of the results presented using our laminar model to the turbulent regime requires careful consideration. First, it is evident that the magnitude of the heat transfer is expected to increase in the case of turbulence as a turbulent flow has improved heat transfer properties due to its thin boundary layer. However, turbulence might also induce additional circulation zones, located at the stator periphery, or sharp edges in general, as shown by Yuan et al. [24]. If located at the core, these circulation zones will reduce the heat transfer drastically in the same matter as the currently observed circulation reduces the heat transfer at the core. Hence, the heat transfer behavior for a turbulent flow cannot be predicted with sufficient confidence using the current model. Future research on the turbulent heat transfer inside the AFPM machine is recommended.

6. Conclusions

A multi-physics, spline-based, finite element model has been developed and validated to analyze the heat and flow behavior of Axial-Flux Permanent Magnet (AFPM) machines. The employed Isogeometric Analysis (IGA) with local refinement capabilities has shown to yield accurate approximations for the considered problems.

The simulations conducted using this model have provided new insights into the laminar heat transfer of open-end rotor-stator configurations, most prominently:

• When dealing with a partially heated stator surface, different heat transfer behavior is observed compared to the literature based on a uniformly heated stator surface. By partially heating the stator, one enables the entering fluid to dissipate its thermal energy in its path toward the axis of rotation. Decreasing the gap-length enhances this effect, as a lower surface-area-to-volume ratio is obtained, thereby effectively increasing the dissipation rate.
• Adding air-inlets will improve the overall cooling of the stator, but will locally worsen it inside the gap near the core surface. The main reason for this is the circulation zone at the core surface, which increases in size with an increase in gap-length or Reynolds number. It is noted that this observation for the laminar flow case does not necessarily translate to the turbulent regime.

When cooling is to be optimized inside an AFPM, the following should be considered:

• An air-inlet is always desired, even though the heat transfer at the core may locally become worse. The cooling through the remaining sides of the stator is expected to be substantially larger compared to having no air-inlet.
• A superposed throughflow at the air-inlet is desired over a conventional open air-inlet configuration. The reason for this is that one can then influence the flow field and therefore control the temperature of the core to some extent. The superposed flow is also able to remove the troublesome circulation zone near the core, increasing the heat transfer.
• An additional cooling circuit that could actively circulate a liquid, i.e., water, refrigerants, etc., will benefit the overall cooling of the stator. The circuit should be positioned such that it can dissipate the heat from the core both efficiently and practically, i.e., preferably located at the back of the stator disk.
• One could use highly conductive non-metal materials for the stator and rotor, dissipating the heat from critical areas even faster.