Flow Resistance Modeling for Coolant Distribution within Canned Motor Cooling Loops

Taylor–Couette–Poiseuille (TCP) flow dominates the inner water-cooling circulation of canned motor reactor coolant pumps. Current research on TCP flow focuses on torque behaviors and flow regime transitions through experiments and simulations. However, research on axial flow resistance in a large Reynolds number turbulent state is not sufficient, especially for the various flow patterns. This study is devoted to investigating the influence of annular flow on the axial flow resistance of liquid in the coaxial cylinders of the stator and rotor in canned motor reactor coolant pumps, and predicting the coolant flow distribution between the upper coil cooling loop and lower bearing lubricating loop for safe operation. The axial flow resistance, coupled with the annular rotation, is experimentally investigated at a flow rate with an axial Reynolds number, Rea, from 2.6 × 103 to 6.0 × 103 and rotational Reynolds number, Ret, from 1.6 × 104 to 4.0 × 104. It is revealed that the axial flow frictional coefficient varies against the axial flow rate in linear relation sets with logarithmic coordinates, which shift up when the flow has a higher Ret. Further examination of the axial flow resistance, with the Rea extending to 3.5 × 105 and Ret up to 1.6 × 105, by simulation shows gentle variation rates in the axial flow frictional coefficients against the Rea. The relation curves with different Ret values converge when the Rea exceeds 3.5 × 105. A prediction model for TCP flow consisting of a polygonal approximation with logarithmic coordinates is developed to estimate the axial flow resistance against different axial and rotational Reynolds numbers for the evaluation of heat and mass transfer during transition states and the engineering design of the canned motor chamber structure.


Introduction
Canned motor pumps have been used in third-generation nuclear power plants, such as the Westinghouse AP1000 advanced passive plant [1] and the SNPTC (State Nuclear Power Technology Corporation) CAP1400 advanced passive plant [2], to circulate primary reactor coolant throughout the reactor core. The level of security has been greatly improved by the elimination of the seal leakage risk brought about by changing the dynamic seal to a static seal when the pumping liquid is introduced into the inner clearance of the canned motor. The fluid in the gap forms an internal cooling system in the canned motor and circulates under the driving of the auxiliary impeller [3], as shown in Figure 1.
The fluid in the motor clearance plays an important role in the cooling of windings and lubrication of bearings. The coolant flow distribution between the upper coil cooling loop and lower bearing lubricating loop depends on the flow resistance.
Since the area of the motor can and flywheel take up the main part of the pump, the dominant pattern of the inner coolant flow is Taylor-Couette-Poiseuille (TCP) flow, which is described as the axial liquid flow between two concentric cylinders with a rotating inner surface, as (2020) 33:14 shown in the left portion of Figure 2. TCP flow is typically used for studying the stability theory of fluid dynamics, momentum transport, boundary layer heat transfer behavior, and the transition of flow regimes. The axial resistance of TCP flow is another research topic that is important in the design of canned motors since the annular gap decides the flow distribution in the upwards and downwards channels. Yamada made a theoretical formula to estimate the resistance of flow through an annulus with an inner rotating cylinder considering the assumed velocity distribution in the gap and validated it experimentally under limited hydraulic parametric conditions [4]. Nouri investigated the laminar and turbulent flow of Newtonian and non-Newtonian fluids in a concentric annulus with rotation of the inner cylinder [5]. He found that, when the dimensionless Rossby numbers were similar, the swirl velocity profiles were the same, and, therefore, the axial drag coefficients were also similar. Kim measured the pressure losses of different working fluids in TCP flow and drew the conclusion that the increase in flow disturbance caused by the Taylor vortex resulted in an increase in the skin friction coefficient, but he did not provide an in-depth explanation [6][7][8]. Kristiawan investigated the components of the wall shear rate via experiments and found that the axial distributions of the wall shear rate components averaged over the perimeter were similar to the distribution in steady Taylor vortices [9]. However, Kristiawan's tests were performed in slow axial flow, which cannot be introduced in the turbulent state. Huisman studied the velocity distribution in TCP flow by laser doppler anemometry and corrected the curvature effect via the use of a ray-tracer [10]. However, his research subject was flow with an outer cylinder rotation. Hashemian predicated the velocity profiles and frictional pressure losses in annular yield-power-law (YPL) fluid flow by numerical methods, and found that the effect of eccentricity was more significant to the reduction of the pressure loss than the radius ratio [11]. Aubert also investigated TCP flow experimentally through laser doppler velocity (LDV) measurements [12]. He put the emphasis on the velocity profiles and heat transfer characteristics, but failed to discuss the axial pressure losses. Yew Chuan measured the pressure loss in the rotor-stator gap [13] and presented the entry effect of air flow entering the annular gap with inner cylinder rotation [14]. Yew's study demonstrated that the flow resistance caused by rotation could be significant, and concluded that the friction factor of laminar flow was not affected by rotation. Sun Chao summarized measuring techniques for turbulent Taylor-Couette flow and provided an experimental point of view on high precision experimental setups [15]. However, the experiments in the review had no axial flow, which indicated that the TCP flow measurements were far from comprehensive and mature. Nouri-Borujerdi experimentally studied the friction factor and heat transfer behavior of TCP flow with a smooth and slotted surface. The results showed that the slot depth enhanced the friction factor, especially at higher effective Reynolds numbers [16], and optimal geometric parameters for the channel were recommended [17].
Most of the studies on TCP flow mentioned above used experimental methods, while numerical methods have been adopted with the progress of computer science in recent years. As early as 1998, Azouz made an  (2020) 33:14 evaluation of three turbulence models for turbulent flow in concentric annuli [18]. He observed that the one-layer mixing length model, two-layer mixing length model, and two-equation model performed similarly when simulating TCP flow. With the improvement of computer performance, direct numerical simulations (DNS) and large eddy simulations (LES) have become advanced methods for studying fluid dynamics with high precision. Rodolfo Ostilla-Mónico is a representative scholar that has performed DNS on Taylor Couette flow from the aspects of angular momentum transfer [19], boundary layer dynamics [20], radius ratio considerations [21], and the effects of domain size [22]. In addition, he provided flow features in detail. Rodolfo also explored the large-scale structure of Taylor-Couette turbulence through LES and regarded it as a useful tool for fast exploration to check for the presence of axially pinned large-scale structures [23]. Akihiro Ohsawa studied the through-flow effects of TCP flow with Reynolds numbers that varied from 500 to 8000 through LES, and found that the friction factor in the axial direction varied with the flow state [24]. Dhaval Paghdar investigated the effects of angular velocity and eccentricity on Taylor-Couette flow and revealed the conditions for the occurrence of Taylor vortices [25]. Though DNS and LES provide more accurate solutions in the time domain, their limitation is obvious; that is, the long calculation time required when the Reynolds number is large leads to less application in engineering design. Compared with the DNS and LES, the Reynolds Average Navier-Stokes (RANS) method has advantages in computing efficiency and a wider range of dynamic parameters. Jacobs compared the results of RANS models with DNS and experiments and found that all models, except for the standard k-ε model, predicted mass flow rates that agreed with the experimental values, but failed to accurately predict the near-wall turbulence dissipation rate [26]. To get more accuracy from the shear stress transport (SST) k-ω model, Dhakal made a modified model, which produced superior performance, especially in the presence of rapid rotation or strong streamline curvature [27,28]. Neto also compared several RANS models in TCP flow research and indicated that the results of the RANS models were very close to each other. However, the SST k-ω model and the Reynolds stress model (RSM) showed slightly better predictions in some aspects [29]. David Shina modeled wide gap Taylor-Couette flow by using an implicit finite volume RANS scheme with a realizable -model [30] and validated the numerical method by experiments [31,32]. Through the above-mentioned studies, RANS models were proved appropriate for simulating TCP flow in a highly turbulent state. Among the RANS models, the SST k-ω and RSM models are the most accurate.
As previously mentioned, numerous studies on the Taylor-Couette system without axial flow between the two coaxial cylinders have been performed experimentally and numerically. However, when additional axial flow was present, existing discussions became inadequate, especially in highly turbulent states in both the axial and circular directions. The rotational speed of the Reactor Coolant Pump (RCP) in the CAP1400 was approximately 1450 r/min, and the peripheral velocity of the inner can was high due to the large inner radius of the annulus, which made the rotational Reynolds number reach up to 5.0 × 10 5 . Furthermore, the mass flow rate of the water in the narrow canned motor gap was above 35 kg/s, which made the axial Reynolds number reach up to 4.2 × 10 4 . The coupled turbulent flow puts forward a challenge in predicting the axial flow resistance of TCP flow.
In this study, the axial flow resistance characteristics of TCP flow in a large Reynolds number turbulent state in a canned motor was investigated via numerical simulation and experiments. A periodic calculation domain and SST k-ω turbulence model was applied in the numerical modeling, which was validated by Yamada's experiments [4], Aubert's measurements [12], and the designed experimental test. Using simulations and experiments, a simplified model was proposed to calculate the axial frictional coefficient accompanied by TCP flow. The investigation provides the basis for the inner cooling clearance design of the canned motor pump and helps to avoid unstable flow to ensure steady operation of the motor.

Geometry, Flow Governing Equations, and Boundary Conditions
In this study, the typical TCP flow style is considered as shown in Figure 2, in which the inner cylinder is the only rotating wall. The geometry of TCP flow is characterized by three parameters: the radius of the inner cylinder, r i , the radius of the outer cylinder, r o , and the axial length, L . The inner wall has angular velocity, ω , and the flow in the gap can be described by the following parameters: The characteristic dimension, d = r o − r i , the radius ratio, η = r i r o , the axial Reynolds number, Re a =Ū a d ν , and the rotational Reynolds number, Re t = r i ωd ν , where ν denotes the kinematic viscosity of the fluid, and Ū a represents the mean axial velocity.
The fluid in the annular gap was water, and the thermal properties of the water were fixed at 30 °C, which was the same as the experimental tests. The flow was considered to be steady-state and incompressible in the isothermal In the equations, ρ represents the fluid density, P is the pressure, U is the velocity, and f is the source term.
According to the work of Rodolfo Ostilla-Mónico [22], the computational domain could be simplified by using a rotational periodic boundary condition, as shown in Figure 2, to reduce the calculation amount. The assumption of periodicity implied that the velocity components repeated themselves in either the axial or rotational direction. When considering the rotational direction, the pressure was also periodic, but in the axial direction, the pressure was not periodic. Instead, in the axial direction, the pressure drop between modules was periodic. The treatments of the axial pressure drop refers to the method offered by Ansys [33]. In this study, the angle, α , between two periodic sections was set to 15°. The inner wall was designated as the rotating wall, and the outer wall was designated as the stationary wall.

Turbulence Model and Solver Definitions
An efficient way to solve turbulence issues is to use the turbulent model. The most widely used turbulent models are generated from the Reynolds Average Navier-Stokes (RANS) equations, which introduce a Reynolds stress term as shown below: Eq. (3) and Eq. (4) are called Reynolds Averaged Navier-Stokes (RANS) equations and have the same form as the instantaneous Navier-Stokes equations, but the velocities and other variables represent timeaveraged values. The last term of Eq. (4) represents the effects of turbulence and the term ρu ′ iū ′ j must be modeled to close the equation. Two-equation turbulence models have the advantages of robustness, economy, and reasonable accuracy for a wide range of turbulent flows, which makes them popular in industrial flow and heat transfer simulations.
The SST k-ω model in Ansys-Fluent, which was proposed by Menter [34], was modified and applied to close the equation and simulate the flow dynamics of the TCP system. A limiter to the eddy-viscosity formulation was set in the model to obtain proper transport behavior. This limiter can be used to avoid over-prediction of the eddy-viscosity [26]. This feature also made the SST kω model more accurate and reliable for a wider class of flows, such as adverse pressure gradient flows, airfoils, and transonic shock waves. The works of Viera Neto, Jacobs, and Dhakal [26][27][28][29] mentioned in the introduction also showed the rationality of using the SST k-ω model.
A coupled scheme of pressure-velocity coupling was chosen for the study, and significantly improved the rate of solution convergence. PRESTO was adopted as the pressure scheme in the spatial discretization, and is an alternative second order scheme that is often useful when strong body forces exist. A third-order convection scheme conceived from the original MUSCL (Monotone Upstream-Centered Schemes for Conservation Laws) [35] was applied to the momentum, turbulent kinetic energy, and turbulent dissipation rate discretization. The least squares cell-based scheme was used in the gradient spatial discretization, which was less expensive and selected as the default gradient method in the Fluent solver.

Meshing and Grid Independent Analysis
The computational domain was created based on the experimental facility and meshed with the ICEM CFD software. The radii of the inner and outer cylinders were 76.5 mm and 79 mm, respectively. The gap thickness, d, was 2.5 mm, and the axial length, L, was 25 mm (i.e., ten times the d). The rotational speed of the inner cylinder was set to 1231.2 r/min, and the mass flow rate was 0.1 kg/s per circumferential period. The temperature of the water was set to 30 °C, the density, ρ, was 995.37 kg/m 3 , and the dynamic viscosity, μ, was 7.49 × 10 −6 Pa·s.
In the simulation, three main parameters were used to describe the grid density, namely the node numbers in the radial, circumferential, and axial directions, as shown in Figure 3. The value of y+ near the boundary wall was set to a value less than three to meet the requirements of the turbulent model. The grid inflation ratio in the radial direction was set to 1.1. Four sets of grids were designed to analyze the sensitivity. Table 1 shows the mesh parameters and iteration times. The calculations were performed on a computer with an Intel i7-6700k CPU, and 32 GB of RAM.
The axial shear stress, τ, and the axial velocity distribution were chosen for the grid independence analysis. The axial shear stress, τ, is defined by Eq. (5): In Eq. (5), F a is the axial force induced by the water flow at the wall, A is the surface area, and i/o refers to the inner or outer cylinder.
The non-dimensional axial velocity is defined by Eq. (6): In Eq. (6), U a is the physical axial velocity, and U τ is the axial wall shear rate defined by Eq. (7): The non-dimensional wall distance, r + , is defined by Eq. (8): In Eq. (8), r is the radial coordinate, and r i/o is the radius of the inner or outer cylinder.
The axial shear stresses of the four mesh sets in Table 1 were calculated with the results shown in Figure 4(a). It was indicated that the shear stress increased with the grid density until it reached a steady state. The axial velocity distributions along the radial gap are represented in Figure 4(b). From the partial enlargement of the curves, the results of mesh 3 and mesh 4 with high mesh density varied slightly. Considering both computational efficiency and accuracy, the third mesh density was adopted to simulate the flow.

Experimental Regime and Method
In order to identify the axial flow resistance in the gap flow of two concentric cylinders with a rotating inner wall, an inner flow test rig for canned motors was designed as shown in Figure 5. The experimental device consisted of two main parts, a rotor shaft and a shell, both made of stainless steel. The shaft was driven by a motor, and had a rotational speed ranging from 500 r/ min to 1500 r/min. A tachometer was employed between the shaft and rotor to measure the motor revolutions. Ten pressure sensors were placed along the axial direction to measure the axial pressure drop. A pump was used to fill the equipment with water, adjust the flow rate, and maintain a relatively high pressure to prevent cavitation. A   (2020) 33:14 flowmeter installed in the pump export pipeline was used to measure the volumetric flow rate. The axial flow frictional coefficient is defined by Eq. (9): where Ū a is the average axial velocity calculated by Q V S , in which Q V is the volumetric flow rate and S is the cross-sectional area of the annulus. The axial distribution of static pressure on the outer cylinder was shown to be linear, which was in accordance with the work of Nouri [36]. The ratio dP dL , which is the pressure loss along the axial direction, can be determined by Eq. (10) via experiments.
In Eq. (10), the measured pressure values of sensor 2, sensor 3, and sensor 4 were used to mitigate the unsteady values of sensors 1 and 5 due to end effects. ΔL is the space between sensors and g is the gravitational acceleration.
Based on the experimental data of the flow rate, the pressure (in terms of Eqs. (9) and (10)), and the relationship between the axial flow resistance C f and the axial flow rate Re a can be obtained.

Experimental Results and Validation of the Simulation Method
The axial pressure-drop performance of the TCP flow was tested at five rotational speeds, which were transformed into rotational Reynolds numbers. The axial flow rates were set in the range of 0.96-2.16 kg/s, and transformed into axial Reynolds numbers. Both the experiments and simulations were carried out under the same parameters.
The results of the designed experimental test and simulation are shown in Figure 6, where the scatter points represent the experimental data and the dashed lines represent the computed results. With a fixed axial flow rate, the frictional coefficient increased with the rotational Reynolds number, which indicated that the cylinder rotation could lead to an increase in the axial flow resistance. It was also observed that the flow frictional coefficient decreased against the axial Reynolds number, and the decline became steeper with the increase in rotational Reynolds number.
The comparison displayed in Figure 6 showed good consistency between the experiments and simulations. Due to the test rig limitations, the radius ratio and axial Reynolds number were not large enough to make an analogy to the working situation of the RCP. In order to validate the numerical method for a wider scope, the experimental results of Yamada [4], with different radius ratios of 0.897 and 0.971 and a larger axial Reynolds number, were further compared with the results of the simulation. The comparison between the simulation and Yamada's test is given in Figure 7. The simulated results agreed well with the experiments for different radius ratios and a larger axial Reynolds number, which indicates that the numerical method was effective in predicting the axial frictional coefficients of TCP flow.
In addition to macro validation through the curves representing the relationship between C f and Re a , the simulated velocity profile was compared with the test results from Aubert [12] by means of a two-component LDV measurement. The normalized axial velocity and radius defined by Eq. (11) and Eq. (12) were used in the comparison: (2020) 33:14 The comparison of the axial velocity distribution from the simulations and Aubert's experiments are given in Figure 8. As shown in the figure, the flow divided into three regions: two boundary layers and a central region, and the simulated data were in good agreement with the measurements in the central region. The test was invalid in the boundary layers due to device constraints such as the wall curvature.
A conclusion could be drawn from the above analysis that the numerical method, with periodic boundary conditions and the SST k-ω turbulence model, was effective in modeling and simulating TCP flow over a wide range when the rotational Reynolds number was in the turbulent state.

Factors Affecting the Axial Flow Resistance in TCP Flow
Two parameters, namely, the radius ratio and rotational Reynolds number, noted in the test and calculation, were investigated in detail by numerical methods.
Firstly, the radius ratio effect was examined. A total of four annuluses were designed based on the dimensions of the experimental device, with the same inner cylinder radius but various outer cylinder radii. The parameters are illustrated in Table 2.
The frictional coefficients under two rotational speeds were calculated in each case, as shown in Figure 9. It can be seen that the lines clearly gathered into two groups with different rotational Reynolds numbers, symbolizing the influence of the rotational speed. With the increase in axial flow rate, the frictional coefficient decreased steeply and tended to be steady. Within the radius ratio range from 0.950 to 0.987, the influence of the radius ratio on the frictional coefficient, C f , could be ignored.
Secondly, the effect of cylinder rotation was analyzed. The axial velocity profile was extracted, since the velocity gradient directly impacted the flow resistance. Seven sets of simulations were performed based on the experimental parameters, and the results are presented in Figure 10. The axial velocity in Figure 10(a) was normalized by Eqs. (11) and (12), and Figure 10(b) was nondimensionalized by Eqs. (6) and (8).
Due to the rotation of the inner wall, the flow field underwent great changes compared with the non-rotating situation. It can be seen from Figure 10(a) that the   This resulted in the growth of the velocity gradient near the wall. This change can be explained by the transport theory of Paoletti [37] and Brauckmann [38], where circumferential shear motion brings in extra energy and increases the intensity of the turbulence. The velocity profile, compared with the wall law proposed by Coles [39], is illustrated in Figure 10(b). The axial velocity profile, which was consistent with the wall law when there was no cylinder motion, deviated from the log law with an increase in cylinder rotation. A larger rotational Reynolds number resulted in a bigger deviation due to the additional tangential component. The conclusion can be drawn that the axial flow was hindered by the rotary motion and had to overcome the resistance of the external energy.

Prediction Model of the Axial Resistance in TCP Flow
As shown in Figure 6, Figure 7, and Figure 9, the relationship curves between the axial flow frictional coefficient and axial Reynolds number displayed two noticeable trends.
1) The frictional coefficient decreased and tended to be steady with the increase in axial Reynolds number; 2) The frictional coefficient could be considered infinitely large when the axial Reynolds number was close to zero.
These features make the curve appropriate for fitting by the power function. Suppose the curves meet the description of Eq. (13), in which the parameter B(Re t ) symbolizes a function of the rotational Reynolds number, and C is the constant exponent. To simplify the solution of coefficients C and B, the natural logarithm of both sides of Eq. (13) is taken as denoted by Eq. (14): The experimental data in the rotary situation were also transformed by the natural logarithmic method and are presented in Figure 11. The points could be fitted linearly, and the slopes of the linear fittings appeared to be equal for different rotation speeds. The expression of the fitting function could be derived through the least squares method and is expressed in Eqs. (15) and (16), indicating the influence of the axial Reynolds number and rotational Reynolds number on the frictional coefficient in axial flow.
The established prediction models presented in Eqs. (15) and (16) were developed from the experimental results and subjected to certain axial Reynolds numbers. A total of four sets of simulations were carried out with a wider range of Re a values as shown in Figure 12. The dashed lines in the picture denote the results calculated with the empirical prediction model as expressed in Eqs. (15) and (16), and the lines with differently shaped points represent the simulated results from the Fluent software.
Three distinct trends occurred with the increase in axial Reynolds number, as shown in Figure 12, and included a rotation dependent region, a transition region, and an ultimate region. In the rotation dependent region, cylinder rotation played an important role in affecting the (13) Figure 11 Logarithmic results of the experimental data and fitting flow resistance, while in the ultimate region, the curves converged together, and cylinder rotation had no apparent effect. This partitioning feature was consistent with the modes of the flow regimes that were proposed by Kaye Joseph [40]. Correspondingly, the flow regime in the rotation dependent region was turbulent flow plus Taylor vortices, and, in the ultimate region, was fully developed turbulent flow. In the ultimate turbulent region, the forecast formula denoted by Yamada, was relatively accurate compared with the simulation. By combining the empirical prediction model, as expressed in Eqs. (15) and (16), under the testing condition with the experimental results of Yamada from the higher axial flow rate region, the empirical model (illustrated in dashed lines) could be used to predict the relationship between the frictional coefficient, C f , and the axial Reynolds number, Re a , as well as the rotational Reynolds number, Re t , with slight differences from the simulated relations with the Fluent software. The transition point (Re a ′, C f ′ ) could be determined from Eqs. (16) and (17) From what has been discussed above, the axial flow resistance of TCP flow in the RCP could be calculated by Eqs. (16) or (17). Taking the canned motor of the CAP1400 RCP as an example, when the axial Reynolds number was about 4.2 × 10 4 and the rotational Reynolds number was around 4.4 × 10 5 , the logarithmic value of the frictional coefficient was − 3.8 in terms of Eq. (16). Simulations were also performed to investigate the axial flow resistance of TCP flow in the RCP at rated rotational speeds, and a comparison between the simulation and (17) C f = 0.0675Re −0.24 a , the prediction model is presented in Figure 13. It can be seen from the figure that the calculation results from the two methods were close, and the error percentage was within 6%. However, the prediction from the empirical model was much more efficient.

Conclusions
The axial flow resistance characteristics of TCP flow were investigated by both experimental and numerical methods. The following conclusions were drawn.

Based on experiments analyzing the relation-
ship between the frictional coefficient, C f , and axial Reynolds number, Re a , as well as the rotational Reynolds number, Re t , an empirical prediction model, C f = exp 1.234 · Re 0.135 t · Re −1.027 a , was developed. Combined with the empirical formula C f = 0.0675Re −0.24 a of Yamada with a large axial Reynolds number, an empirical prediction model, in the form of a polygonal approximation with logarithmic coordinates and a wider range of axial Reynolds numbers could be used to achieve a high efficiency prediction with almost the same accuracy as the numerical simulation. 2. The cylinder rotation exhibited a strong influence that increased the axial flow resistance at low axial flow rates, which degraded until invisible with increasing axial Reynolds numbers. 3. The radius ratio had less effect on the axial flow resistance compared with the influence of the rotational Reynolds number. 4. The unstable flow region, as indicated in the transition area of the empirical prediction model, should (2020) 33:14 be prevented in the parameter design of the inner cooling clearance.