Design Optimization of Balancing Holes in a Centrifugal Pump Master’s thesis in Applied Mechanics BALA KUMARESH THILEEP KUMAR DEPARTMENT OF MECHANICS AND MARITIME SCIENCE CHALMERS UNIVERSITY OF TECHNOLOGY Gothenburg, Sweden 2025 www.chalmers.se www.chalmers.se Master’s thesis 2025 Design Optimization of Balancing Holes in a Centrifugal Pump BALA KUMARESH THILEEP KUMAR Department of Mechanics and Maritime Sciences Division of Fluid Dynamics Chalmers University of Technology Gothenburg, Sweden 2025 Design Optimization of Balancing Holes in a Centrifugal Pump BALA KUMARESH THILEEP KUMAR © BALA KUMARESH THILEEP KUMAR, 2025. Examiner: Håkan Nilsson Professor at Division of Fluid Dynamics Department of Mechanics and Maritime Sciences Chalmers University of Technology Industrial Supervisor: Arash Soltani Christian Wollblad Xylem Water Solutions Sweden Master’s Thesis 2025 Department of Mechanics and Maritime Sciences Division of Fluid Dynamics Chalmers University of Technology SE-412 96 Gothenburg Telephone +46 31 772 1000 Cover: Streamline visualization generated in ParaView, illustrating internal jet-like flow behavior through balancing holes in a semi-open centrifugal pump impeller. Streamlines are colored by relative flow angle magnitude (β). Typeset in LATEX, template by Kyriaki Antoniadou-Plytaria Gothenburg, Sweden 2025 iv Design Optimization of Balancing Holes in a Centrifugal Pump BALA KUMARESH THILEEP KUMAR Department of Mechanics and Maritime Chalmers University of Technology Abstract Centrifugal pumps are widely used in fluid transport systems due to their ability to deliver high flow rates in a compact, low-maintenance design. Despite their well- defined performance characteristics, a major operational challenge is the generation of axial force, a force imbalance caused by the pressure difference across the front (suction) and back (pressure) sides of the impeller. If the axial force is not prop- erly investigated during the design process, it can cause bearing failure, increase maintenance requirements, and significantly reduce both hydraulic efficiency and pump lifespan. Balancing holes are commonly introduced to mitigate axial force by reducing the pressure imbalance. While effective in reducing axial force, balancing holes introduce fluid recirculation from the pressure side to the suction side of the impeller, leading to a reduction in hydraulic performance. Therefore, identifying an optimal balancing hole configuration is essential to minimize axial force without compromising overall efficiency. This study investigates the development of a multi-objective target function us- ing both scalarization and Pareto-based methods to identify an optimal balancing hole design configuration that reduces axial force while maintaining the hydraulic performance. CFD simulations were carried out using a steady-state approach with implicit local time stepping, in conjunction with a non-linear k−ε turbulence model and the frozen rotor approach, to capture the internal flow behaviour and to extract performance metrics—head, power, and axial force—for all generated configurations during the optimization process. To explore the design space effectively, Sobol se- quence sampling was used to generate a diverse set of initial configurations. The target function was formulated by normalizing head and power at the Best Effi- ciency Point (BEP) and axial force across part-load, BEP, and over-load conditions, then iteratively optimized using Bayesian optimization and target function scoring to identify an optimal balancing hole design. A Pareto front analysis was finally conducted to analyze the trade-offs between hydraulic performance and axial force reduction for the generated configurations throughout this study. The findings revealed that the formulated target function using scalarized method effectively guided the optimization process toward configurations that reduced axial force by up to 40.59 % while maintaining head and power within 99.98 % and 99.87 % of the reference, respectively. Iterative Bayesian optimization and Pareto front analysis consistently led to a key design feature, which was also supported by the flow field visualization. Keywords: Balancing hole, Multi-objective target function, Sobol sequence sam- pling, Bayesian optimization, Pareto front analysis. v Acknowledgements First and foremost, I would like to express my heartfelt gratitude to my supervisors, Arash Soltani and Christian Wollblad, at Xylem for their continuous support, en- couragement, and guidance throughout my master’s thesis. From helping me with fundamental details to complex decision-making, their support has been truly in- strumental in every stage of this journey. Their mentorship went beyond technical supervision—they were always available and patient, and I am deeply thankful for the knowledge and experience I gained working alongside them. I would also like to sincerely thank Fredrik Söderlund, my manager at Xylem, for providing me this incredible opportunity and for fostering a supportive environ- ment throughout my time at the company. I am grateful to Tobias Pilstrand for his valuable assistance with Bayesian optimization techniques and for always being approachable when I needed help. A special thanks to Viktor Bredvad for providing the impeller geometry and helping me understand the design aspects involved—his support enabled me to create multiple geometries in a short time, which was crucial for generating a wide range of configurations used in this thesis. My thanks extend to the entire Xylem team for welcoming me and for their tech- nical insights, discussions, and encouragement, which made my thesis work both meaningful and enjoyable. I would also like to thank Chalmers University of Technology for providing an ex- cellent academic environment and resources that supported my learning and growth throughout the program. I am grateful to Professor Håkan Nilsson, my examiner, for his valuable feedback and for supporting the academic evaluation of this thesis. A special mention goes to Professor Hudong Yao, whose referral played a key role in helping me secure this thesis opportunity. I would also like to thank the entire De- partment of Mechanics and Maritime Sciences and all the professors whose courses helped build the foundation for this work. Finally, I wish to thank my family and friends for their constant encouragement, support, and belief in me during this wonderful journey. And last but not least, I would like to acknowledge myself—for staying committed, believing in my abilities, and pushing through every challenge along the way. Bala Kumaresh Thileep Kumar, Gothenburg, June 2025 vii List of Acronyms Below is the list of acronyms that have been used throughout this thesis listed in alphabetical order: BEP Best Efficiency Point BH Balancing Hole CDF Cumulative Distribution Function CFD Computational Fluid Dynamics CRPS Continuous Ranked Probability Score EI Expected Improvement PDF Probability Density Function RSE Relative Squared Error TF Target Function ix Nomenclature Below is the nomenclature of symbols that have been used throughout this thesis. Symbols Centrifugal Pump Ẇc Impeller power (work rate on the fluid) τA Torque acting on the impeller Ω Angular velocity of the impeller ṁ Mass flow rate U1 Blade speed at inlet U2 Blade speed at outlet cθ1 Tangential component of the absolute velocity at inlet cθ2 Tangential component of the absolute velocity at outlet ∆W Specific work (work done on fluid per unit mass) Htheo Theoretical head developed by the pump H Total head developed by the pump g Acceleration due to gravity PShaft Shaft power PHydraulic Hydraulic power imparted to the fluid ρ Density of the fluid η Hydraulic efficiency of the pump Q Volumetric flow rate Q∗ Normalized flow rate QBEP Flow rate at BEP condition P2 Pressure at outer radius of the impeller Hp2 Static pressure C2 Absolute velocity of the liquid at impeller outlet xi F Axial force acting on the impeller T Net axial force acting along the pump shaft F1 Pressure force on the rear shroud F2 Pressure force on the impeller hub F3 Force acting across the blade tip clearance F4 Dynamic reaction force at the impeller inlet AClearance Clearance area of the impeller DMean Mean diameter of the impeller s Radial clearance between rotating and stationary components Dmin Minimum diameter of the balancing hole Dmax Maximum diameter of the balancing hole ATotal Total area of the balancing holes Aper, hole Area of balancing hole Computational Fluid Dynamics k Turbulent kinetic energy ε Turbulent dissipation rate µ Dynamic viscosity vi Instantaneous velocity component p Instantaneous pressure νt Eddy viscosity Sij Mean strain rate tensor aij Anisotropic part of the Reynolds stress tensor δij Kronecker delta Multi-Objective Optimization N Number of samples fi(xj) Value of ith objective for jth sample fi,ref Reference value for the ith objective ∆i Mean absolute deviation Ii Inverse of the mean absolute deviation wi Normalized weight assigned to ith objective m Number of objectives xii Ik Inverse deviation of kth objective fi(x) ith objective function F (x) Scalarized objective function yref, i Reference value for ith data point ŷi Observed or simulated value for the ith data point yref Mean of all reference values n Total number of data points yi Observed or simulated value for the ith configuration f(x) Objective function µ(x) Mean function E[f(x)] Expected value K(x, y) Covariance function (or kernel) k Variance parameter σ Length-scale parameter of the function Z Standardized improvement f ∗ Best value observed Φ(Z) Cumulative distribution function ϕ(Z) Probability density function DN Discrepancy of the first N Sobol points θ Model parameters D Observed dataset x∗ New input y∗ Predicted output at the new input x∗ Xi Simulated result of the quantity (Head or Power) for ith impeller configuration with balancing hole at BEP (Q∗ = 1.0) condition Xref Simulated result of the quantity (Head or Power) for the reference impeller without balancing holes at BEP (Q∗ = 1.0) condition σX Standard deviation of the simulated results (Head or Power) across the 12 configurations generated using Sobol sequence at BEP (Q∗ = 1.0) condition Fi Simulated axial force for the ith impeller configuration with balanc- ing hole at BEP (Q∗ = 1.0) condition Fref, i Simulated axial force for the reference impeller without balancing hole at BEP (Q∗ = 1.0) condition Fi, 0.5 Simulated axial force for the ith impeller configuration with balanc- ing hole at part-load (Q∗ = 0.5) condition xiii Fi, 1.0 Simulated axial force for the ith impeller configuration with balanc- ing hole at part-load (Q∗ = 1.0) condition Fref, 0.5 Simulated axial force for the reference impeller without balancing hole at BEP (Q∗ = 0.5) condition NormalizedP Normalized value of power NormalizedH Normalized value of head NormalizedAF Normalized value of axial force InvP Inverse deviation for power InvH Inverse deviation for head InvAF Inverse deviation for axial force wP Weightage value for power wH Weightage value for head wAF Weightage value for axial force TFPredicted Predicted target function score obtained using Bayesian optimiza- tion TFScalarized Target function score obtained using the scalarized form (∆η)max Maximum efficiency loss ηRef Efficiency for the reference impeller ηi Efficiency for the ith impeller configuration with balancing hole xiv Contents List of Acronyms ix Nomenclature xi 1 Introduction 1 1.1 Problem Statement . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.2 Objective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 1.3 Limitations of the Study . . . . . . . . . . . . . . . . . . . . . . . . . 3 2 Theory 5 2.1 Centrifugal Pump . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.1.1 Working Principle . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1.2 Classification of Centrifugal Pump considered in this study . . 6 2.1.3 Energy Transfer within a Centrifugal Pump . . . . . . . . . . 7 2.1.3.1 Head Developed . . . . . . . . . . . . . . . . . . . . 7 2.1.3.2 Power Consumption & Efficiency . . . . . . . . . . . 8 2.1.4 Performance Characteristics of a Centrifugal Pump . . . . . . 8 2.1.4.1 Q∗ − H Curve . . . . . . . . . . . . . . . . . . . . . 9 2.1.4.2 Q∗ − P Curve . . . . . . . . . . . . . . . . . . . . . . 9 2.1.4.3 Q∗-η Curve . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 Forces Acting on a Centrifugal Pump . . . . . . . . . . . . . . . . . . 10 2.2.1 Axial Force on Semi–Open Impeller . . . . . . . . . . . . . . . 10 2.2.2 Balancing Holes . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.2.2.1 Effect of Balancing Holes on Pressure Distribution and Axial Force . . . . . . . . . . . . . . . . . . . . . 13 2.2.2.2 Design Configurations . . . . . . . . . . . . . . . . . 13 2.2.3 Clearance Area . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.3 Computational Fluid Dynamics . . . . . . . . . . . . . . . . . . . . . 14 2.3.1 Turbulence Modeling . . . . . . . . . . . . . . . . . . . . . . . 16 2.3.2 Standard k − ε Turbulence Model . . . . . . . . . . . . . . . . 16 2.3.3 Cubic k − ε Turbulence Model . . . . . . . . . . . . . . . . . . 17 2.3.4 Multiple Reference Frame model - Frozen Rotor Approach . . 18 2.4 Multi-Objective Optimization . . . . . . . . . . . . . . . . . . . . . . 19 2.4.1 Pareto Optimality . . . . . . . . . . . . . . . . . . . . . . . . 19 2.4.1.1 Preferences and Utility Function . . . . . . . . . . . 20 2.4.2 Normalization . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 xv Contents 2.4.2.1 Relative Squared Error (RSE) . . . . . . . . . . . . . 22 2.4.2.2 Pointwise Relative Squared Error . . . . . . . . . . . 22 2.4.3 Bayesian Optimization . . . . . . . . . . . . . . . . . . . . . . 23 2.4.3.1 A Bayesian Approach to Optimization . . . . . . . . 23 2.4.3.2 Gaussian Process . . . . . . . . . . . . . . . . . . . . 23 2.4.3.3 Acquisition Function . . . . . . . . . . . . . . . . . . 24 2.4.3.4 Bayesian Optimization Work Flow . . . . . . . . . . 25 2.4.4 Sobol Sequence Sampling . . . . . . . . . . . . . . . . . . . . . 26 2.4.4.1 Discrepancy and Uniformity Properties . . . . . . . . 26 2.4.5 Posterior Predictive Distribution . . . . . . . . . . . . . . . . 26 3 Methodology 27 3.1 Design configurations for balancing hole optimization . . . . . . . . . 27 3.1.1 Number of balancing holes . . . . . . . . . . . . . . . . . . . . 28 3.1.2 Balancing hole diameter . . . . . . . . . . . . . . . . . . . . . 28 3.1.3 Circumferential position of balancing hole . . . . . . . . . . . 29 3.1.3.1 Back Side of the Impeller . . . . . . . . . . . . . . . 29 3.1.4 Radial Position of Balancing Hole . . . . . . . . . . . . . . . . 29 3.1.5 Finalization of Design Parameter Ranges . . . . . . . . . . . . 30 3.2 Computational Modeling . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.2.1 Computational Domain . . . . . . . . . . . . . . . . . . . . . . 31 3.2.2 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . 32 3.2.3 Computational Mesh . . . . . . . . . . . . . . . . . . . . . . . 34 3.2.3.1 Baseline Mesh Configuration (Without Balancing Hole) 34 3.2.3.2 Mesh Refinement Study around Balancing Hole . . . 36 3.2.4 Turbulence Modeling and Solver Settings . . . . . . . . . . . . 37 3.3 Initial Sampling for the Optimization Process using Sobol Sequence . 38 3.4 Target Function Development Strategy . . . . . . . . . . . . . . . . . 40 3.4.1 Normalization Strategy . . . . . . . . . . . . . . . . . . . . . . 41 3.4.1.1 Normalization of Head and Power . . . . . . . . . . . 41 3.4.1.2 Normalization of Axial Force . . . . . . . . . . . . . 42 3.4.2 Target Function Formulation . . . . . . . . . . . . . . . . . . . 43 3.4.2.1 Geometric Mean Method . . . . . . . . . . . . . . . 43 3.4.2.2 Scalarized Form with Inverse Deviation Weightage Method . . . . . . . . . . . . . . . . . . . . . . . . . 44 3.4.3 Bayesian Optimization & Pareto Front Analysis . . . . . . . . 45 4 Results 47 4.1 Performance of Reference Impeller (without Balancing Hole) . . . . . 47 4.2 Results of Mesh Refinement Study Around the Balancing Hole . . . . 48 4.3 Performance Evaluation of Sobol-Seeded Configurations . . . . . . . . 51 4.4 Normalization Results of Head & Power for Target Function . . . . . 53 4.5 Normalization Results of Axial Force for Target Function . . . . . . . 53 4.5.1 Baseline Normalization using BEP Reference (Q∗ = 1.0): . . . 54 4.5.2 Worst-Case Normalization across Operating conditions . . . . 54 4.6 Optimization Using the Scalarized Form of the Target Function . . . 55 xvi Contents 4.6.1 First Optimization Cycle using Baseline Normalization of Ax- ial Force . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 4.6.2 Second and Third Optimization Cycle (Worst Case Normal- ization of Axial Force) . . . . . . . . . . . . . . . . . . . . . . 57 4.7 Pareto Front Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.8 Performance Results of Impeller Configurations . . . . . . . . . . . . 62 4.9 Influence of Design Parameters during Optimization Process . . . . . 68 4.10 Comparison of predicted and scalarized values of Target Function . . 69 4.11 Flow Visualization and Interpretation of Balancing Hole Effects . . . 70 4.11.1 Pressure Distribution and Mass Flow Rate Analysis through Balancing Holes . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.11.1.1 BEP Condition — Q∗ = 1.0 . . . . . . . . . . . . . . 71 4.11.2 Jet Interaction and Flow Redistribution through Balancing Hole 73 4.11.2.1 BEP Condition — Q∗ = 1.0 . . . . . . . . . . . . . . 74 5 Conclusion 77 5.1 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78 Bibliography 79 A Appendix I A.1 Python Code for Sobol Sequence Sampling . . . . . . . . . . . . . . . I A.2 Python Code for Pareto Front Analysis . . . . . . . . . . . . . . . . . II xvii Contents xviii 1 Introduction In the modern world of fluid transport, few machines are as indispensable as pumps. They serve as the driving force behind countless industrial and domestic systems, ensuring the efficient and reliable liquid transport through complex networks of pipelines and processing units. From delivering clean water to households and wastewater management to supporting the rigorous demands of industrial opera- tions, pumps play a significant role in contemporary infrastructure. Pumps function by transforming mechanical energy into fluid energy and operate primarily through two mechanisms: fixed-volume displacement (positive displacement pumps) or con- tinuous kinetic energy transfer (dynamic pumps). Among the various types of dy- namic pumps, centrifugal pumps stand out due to their superior performance, com- pact design, and ability to handle high flow rates with relatively low maintenance requirements. Their broad range of applications—from municipal water systems to high-demand industrial processes—underscores their critical role in modern engi- neering and fluid transport systems. As previously noted, centrifugal pumps operate by converting mechanical energy into hydraulic energy. In this process, the fluid gains kinetic energy through the rotating impeller, which is then gradually transformed into pressure energy as the fluid moves through the pump. Despite their well-defined performance character- istics (like efficiency and handling high flow rates), centrifugal pumps face certain operational challenges, with this study focusing specifically on axial loads. These loads result from a pressure imbalance between the front and back sides of the im- peller. If not addressed properly, axial loads can cause bearing failures, increased maintenance requirements, and a decline in both pump performance and service life [1]. Among the various options available to reduce the axial load, Shepherd [2] consid- ered the use of balancing holes to be the simplest and cheapest option. These holes connect the high-pressure region behind the impeller to the low-pressure region at the front, allowing a portion of the fluid to recirculate [3]. While this pressure equal- ization helps reduce axial force, it comes at the cost of hydraulic performance. The recirculated fluid introduces energy losses, which can negatively affect head genera- tion, power consumption, and overall efficiency [4]. The central engineering challenge lies in optimizing the design of balancing holes to effectively reduce axial load without significantly compromising pump performance. Achieving this balance is complex and depends heavily on the configuration of the 1 1. Introduction balancing holes—including their number, diameter [5], and placement—as well as the intricate internal flow behavior within the pump. As industries continue to demand higher-performing pumps, the refinement of balancing hole design remains a key focus in the development of modern centrifugal pump technology. 1.1 Problem Statement Despite their widespread application, there are no robust methods for optimally designing balancing holes that both reduce the axial load and preserve the overall pump performance. In industries, balancing holes are drilled based on guesswork [6], which leads to lower pump performance without significantly reducing the axial force [4, 5]. These limitations highlight the need for a more advanced and reliable op- timization framework —one that leverages CFD techniques and performance-driven algorithms to optimize the hole configuration for improved pump performance. Previous studies by Xylem have adapted similar techniques to address these chal- lenges, providing valuable insights into the flow behavior and design trade-offs. How- ever, they encountered difficulties in defining a comprehensive multi-objective target function capable of simultaneously minimizing axial load and maintaining critical pump performance metrics. This thesis addresses that gap by formulating a cus- tomized multi-objective target function specifically for the Flygt N-centrifugal pump, designed and developed by Xylem, aiming to overcome those limitations. 1.2 Objective Continuing the work initiated by Xylem, this thesis aims to formulate a multi- objective target function that supports the identification of optimal balancing hole design configurations. Instead of delivering a single optimal design for the pump considered in this study, the primary objective is to establish a consistent frame- work using both scalarization and Pareto-based methods to evaluate and compare multiple design configurations. This formulated function is intended to facilitate fu- ture optimization efforts by targeting maximum possible axial force reduction while maintaining the overall hydraulic performance of the pump. To achieve this, de- tailed CFD simulations will be conducted to analyze the hydraulic performance of the generated configurations throughout the optimization process. Based on these, the objectives of this study are classified into primary and secondary objectives: • Primary Objective: Develop a multi-objective target function that inte- grates the axial force and the key performance metrics (Head, Power Consump- tion & Efficiency) at the Best Efficiency Point (BEP) of the pump. Evaluate the effectiveness of the formulated target function by analyzing the generated configurations at BEP but also across a range of operating condition. Iter- atively refine the target function until the selected configuration consistently maintains the pump’s characteristic performance at all investigated operating conditions while minimizing axial force. 2 1. Introduction • Secondary Objective: Analyze CFD simulation results of the generated configurations, to gain insights into internal flow phenomena with and with- out balancing holes. The focus is to provide a physical insight by evaluating pressure distribution and relevant hydrodynamic effects to identify the rela- tionship between design configurations of balancing holes and their influence on axial force and performance metrics of the pump. 1.3 Limitations of the Study This study is subject to the following limitations based on the defined scope, com- putational resources, and time constraints: • The development of the multi-objective target function and the associated analysis results are specific to the geometry and operating range of the cen- trifugal pump with a vaneless diffuser configuration. Generalization to other pump series or types is beyond the scope of this work. • All simulations are carried out under steady-state condition using the Multiple Reference Frame model with the Frozen Rotor approach. Transient simula- tions are not conducted in this study for CFD simulations. • The study exclusively focuses on the reduction of axial forces using optimiza- tion of balancing hole configuration. Radial forces are not considered. • Cavitation is not considered in this study. The simulations assume single phase, incompressible flow. • Temperature effects and thermal-related fluid properties are not included in the simulations. All analyses are performed under isothermal conditions. • The study does not include structural evaluations such as stress analysis, fa- tigue assessment, or deformation of the impeller or surrounding regions near the balancing holes. • A mesh independence study has been carried out in the region around the balancing holes. A full domain mesh independence study for the computa- tional domain was not carried out, as the baseline mesh without balancing holes provided by Xylem was used without modification. • Validation of the optimized target function results was not possible due to the absence of experimental or real-world data. However, the baseline configura- tion without balancing hole is validated against the transient CFD simulation results previously conducted by Xylem which in turn have been validated against measured performance. 3 1. Introduction These limitations are acknowledged to clarify the scope of the study, which focuses on analyzing the hydraulic performance and internal flow behavior associated with balancing hole configurations under defined and controlled conditions. 4 2 Theory This chapter outlines the theoretical concepts relevant to this study and is organized into four main sections: the fundamentals of centrifugal pump, including the energy transfer process, the development of axial forces in a centrifugal pump, an general overview of CFD and turbulence modeling, and the frame-work for multi-objective optimization. 2.1 Centrifugal Pump Centrifugal pumps are hydraulically driven turbomachines designed to transport flu- ids, primarily liquids, by converting mechanical energy into hydraulic energy through the action of centrifugal force [7]. Their main function is to transfer the fluid by increasing its velocity and pressure. Figure 2.1: Schematic representation of a centrifugal pump. Generated using ChatGPT. A typical centrifugal pump consists of an impeller, suction nozzle, discharge nozzle and pump casing [8] as shown in figure 2.1. The impeller comprises a series of blades, increases both the static pressure and velocity of the incoming fluid as it passes through during rotation [9]. The pump casing (either a diffuser or volute type) encloses the impeller and converts the fluid’s kinetic energy into pressure energy. Additionally, the casing is also equipped with suction and discharge sections to guide 5 2. Theory the flow into and out of the pump efficiently. Despite various types of centrifugal pumps, the fundamental operating principle remains consistent, which is discussed in detail in the following section. 2.1.1 Working Principle The fundamental working principle of a centrifugal pump is based on the centrifu- gal force generated by the rotating motion of the impeller. Fluid enters the pumps through the suction nozzle at the center of the impeller (eye of the impeller). As the fluid enters, the rotational motion of the impeller imparts kinetic energy onto the fluid, causing it to accelerate. Due to centrifugal force, which is constantly acting on the fluid particle, the fluid is pushed towards the outer edge of the impeller, thereby increasing its kinetic energy [7]. The curvature and the orientation—especially the arrangement of impeller blades—guide the fluid that is being pushed by the cen- trifugal force smoothly, thereby minimizing energy loss. Simultaneously, due to the geometry of the impeller channel, part of this kinetic energy is converted into pres- sure energy. By the time the fluid reaches the outlet of the impeller, its velocity is high, and its pressure has also increased substantially [3]. Upon leaving the im- peller, the fluid enters the pump casing. The pump casing is designed in such a way to collect the fluid discharged from the impeller and convert part of the remaining fluid’s kinetic energy into pressure energy, typically using a volute type or diffuser type design [3]. In a volute type, the gradually expanding spiral shape of the pump casing reduces the velocity of the fluid leaving the impeller, thereby increasing the fluid pressure towards the discharge. In contrast, the diffuser casing uses a station- ary set of guide vanes surrounding the impeller to achieve the same effect by by gradually decelerating the fluid. According to Bernoulli’s principle, this reduction in velocity results in an increase in static pressure, allowing for efficient pressure recovery before the fluid exits the pump [9]. 2.1.2 Classification of Centrifugal Pump considered in this study Centrifugal pumps are classified based on several factors, including their application, flow direction, impeller design, and the number of stages [8]. This section describes the classification of the centrifugal pump considered in this study. • Flow Type: Radial: The fluid enters the pump axially, along the shaft axis, and is discharged perpendicularly to the direction of the shaft axis. This is a common configuration for pumps handling moderate to high heads. • Impeller Design: Semi-Open Impeller: The impeller blades are attached to a hub on one side, while the other side is exposed to the flow. • Stage Count: Single Stage: A single impeller is mounted on the pump shaft, meaning all energy transfer occurs within one stage. Single-stage pumps are typically used in applications where moderate pressure rise is sufficient. 6 2. Theory • Shaft Orientation: Vertical: The pump shaft is vertically oriented with re- spect to the ground plane. This configuration is commonly used in applications with limited horizontal space or where submersible installation is required. • Casing Type: Vaneless Diffuser: The pump casing features a vaneless diffuser design, characterized by an open annular passage without stationary guide vanes. 2.1.3 Energy Transfer within a Centrifugal Pump The energy transfer within a centrifugal pump can be mathematically described us- ing fundamental principles of fluid mechanics and rotational dynamics. The conser- vation law, particularly the conservation of angular momentum, relates the amount of torque required to impart a change in angular momentum to the fluid, ultimately leading to the fundamental Euler’s equation for turbomachines [9]. This equation governs the energy transfer between the rotating component and the working fluid within a centrifugal pump. Based on this, Euler’s Equation 2.1 quantifies the amount of work done on the fluid by the rotating impeller. The impeller’s work rate (Ẇc) on the fluid is defined as the product of torque and angular velocity: Ẇc = τAΩ = ṁ(u2cθ2 − u1cθ1) (2.1) where, the blade speed u = ω · r, τA is the torque exerted by the impeller, Ω is the angular velocity of the impeller, ṁ is the mass flow rate, u1 and u2 are the blade speeds at the inlet and outlet respectively, and cθ1 and cθ2 are the tangential components of the absolute velocity at the inlet and outlet respectively. Dividing Equation 2.1 by the mass flow rate gives the specific work (or) work done on the fluid per unit mass. Equation 2.2 is referred to as the Euler’s pump equation [9]. ∆W = Ẇ ṁ = τAΩ ṁ = U2cθ2 − U1cθ1 (2.2) 2.1.3.1 Head Developed The specific work expressed in Equation 2.2 can be expressed in terms of head, which represents the energy transferred to the fluid as an equivalent head. By dividing the specific work from Euler’s equation by the gravitational acceleration, the theoretical head developed by the pump can be expressed as: Htheo = ∆W g = U2cθ2 − U1cθ1 g (2.3) Equation 2.3 represents the theoretical estimate of the head developed by the pump based solely on the impeller geometry and the tangential components of the absolute 7 2. Theory velocity at the inlet and outlet, meaning no energy losses due to friction, turbulence, leakage, or other secondary effects [9]. 2.1.3.2 Power Consumption & Efficiency The power consumption of a centrifugal pump refers to the mechanical energy re- quired to drive the pump shaft. It is referred to as shaft power, and it is greater than the hydraulic power due to the presence of various losses within the pump, such as mechanical, volumetric, and hydraulic losses [9]. Equation 2.4 expresses Pshaft the as: Pshaft = Phydraulic η (2.4) where η is the pump’s hydraulic efficiency. The hydraulic power imparted to the fluid, which represents the useful power transmitted through the fluid to overcome the head [9], is defined by Equation 2.5: Phydraulic = ρgQH (2.5) where ρ is the density of the fluid, g is the gravitational acceleration, Q is the volumetric flow rate, and H is the total head developed by the pump. 2.1.4 Performance Characteristics of a Centrifugal Pump In a centrifugal pump, energy transfer is accomplished by changing the state of motion of the liquid as it passes through the rotating impeller. The magnitude of this energy transfer depends upon the geometry of the impeller and the pump’s operating conditions [8]. While the impeller geometry and the rotational speed of the pump remain constant, the fluid velocity varies with changes in flow rate. Pumps often are subjected to operate under varying conditions to meet different system demands, making it essential to evaluate their performance across a range of flow rates. This evaluation is carried out using performance characteristics, which describe how the pump responds to different flow rates in terms of head development, power consumption, efficiency, and Net Positive Suction Head (NPSH) required to avoid cavitation. In this study, the NPSH is excluded from the analysis, as the focus is on optimizing balancing hole design based on internal flow behaviour and its influence on axial force and hydraulic performance. To enable consistent comparison, the flow rate in the performance curves is normalized relative to the flow rate corresponding to the BEP, as defined in Equation 2.6. Each performance characteristics curve is briefly explained in the upcoming subsections to provide a theoretical background of the parameters evaluated in this study. Q∗ = Q QBEP (2.6) 8 2. Theory Q∗ Pe rfo rm an ce C ha ra ct er ist ic s Power Head Efficiency Figure 2.2: Typical performance chart of a centrifugal pump. 2.1.4.1 Q∗ − H Curve The Q∗−H curve represents the amount of head developed by the pump for different flow rates. As illustrated in Figure 2.2, this curve typically exhibits a negative slope, indicating that the head decreases as the flow rate increases. This decline results from increasing internal energy losses within the pump at higher flow rates. Evaluating this curve helps in understanding the pump’s capability to develop the desired head under different operating conditions. 2.1.4.2 Q∗ − P Curve The Q∗ − P curve describes the power consumed by the pump (i.e., shaft power) at different flow rates. As seen in Figure 2.2, this curve follows a positive trend, where the power consumption increases with an increase in flow rate. This increase is due to the fact that a pump requires higher power to transfer greater volume of fluid . Evaluating the Q–P curve is essential for selecting an appropriately sized motor and estimating power requirements under various operating scenarios. 2.1.4.3 Q∗-η Curve The Q∗-η curve represents the proportion of mechanical energy converted into useful hydraulic energy. As shown in Figure 2.2, this curve typically exhibits a peak: efficiency increases with flow rate up to a maximum point—the BEP—and then declines beyond it. Identifying this peak is crucial for determining the pump’s optimal operating point and maximizing overall performance. 9 2. Theory 2.2 Forces Acting on a Centrifugal Pump Centrifugal pumps, like all rotating turbomachines, are subjected to various fluid dynamic forces. These forces arise inherently from the energy transfer process (via angular momentum exchange between the impeller and the working fluid) as well as from internal flow dynamics (including pressure imbalance, momentum changes and unsteady flow effects induced by blade-passage interactions). In a rotating impeller, these forces act in both radial and axial directions and may also include transient components from such interactions. While all of these forces can negatively impact pump reliability and operational lifespan, this study focuses on axial force due to its critical significance in causing bearing failure and reduction in pump life [1]. The following section presents a de- tailed analysis of axial force generation, emphasizing the fluid dynamic mechanisms involved and the associated mathematical formulations. 2.2.1 Axial Force on Semi–Open Impeller Axial force acting on an impeller, regardless of configuration, arises from pressure imbalance between its front (suction) side and the back (pressure) side of the im- peller. A mathematical understanding of this pressure imbalance begins with the head developed by the pump as described in Equation 2.3. This head represents the total energy per unit weight imparted to the fluid, comprising both the static pressure head (Hp2) and dynamic (velocity) head (C2 2/2g), where C2 is the absolute velocity of the liquid at the impeller outlet. Specifically, the pressure P2 at the impeller outlet is equal to Hp2. This pressure acts in all directions, including within the cavity between the impeller hub and pump casing. However, due to the rotational motion of the liquid in this space, the pressure acting on the back side of the hub is not uniform and varies with the radius. For practical purposes, it is often assumed that the liquid within that space rotates with half the angular velocity of the impeller. Using this assumption, the pressure distribution on either side of the semi–open impeller can be derived by applying Bernoulli’s equation in a rotating reference frame. This pressure distribution varies radially inward towards the pump shaft as illustrated in the upper part of Figure 2.3, due to the reduction in the circumferential velocity resulting from centrifugal effects, and can be expressed as: P (R) = P2 − γ(U2 2 − U2) 8g (2.7) where, P2 is the pressure at the outer radius of the impeller R2, U2 = ωR2 is the circumferential speed at the impeller outlet, U = ωR is the circumferential speed at any arbitrary radius R, and γ = ρg is the specific weight of the liquid. This equation explains how the centrifugal effect of the rotating fluid generates a radially varying pressure distribution on both side of the impeller. However, the pressure distribution on the suction side of the impeller follows a more complex profile, typically taking the 10 2. Theory form of a curve, as illustrated in the lower part part of Figure 2.3. The shape of this curve is heavily influenced by the geometry of the impeller blades. While in theory, the pressure distribution could be computed using the usual approach as discussed above, in practice this is highly challenging. This difficulty arises because pressure within the impeller passages is not uniform at a given radius; instead, it varies along the arc between consecutive blades due to complex internal flow phenomena. For practical estimation, it has therefore been suggested that the pressure variation within the impeller passage may linearly vary with radius or proportionally to the square of the radius. Both assumptions have been evaluated, and while neither produces a perfect balance, they have been proven adequate to prevent excessive bearing overload [8]. Pump Shaft Rear Shroud Impeller Vanes Figure 2.3: Schematic diagram of pressure distribution on a semi-open impeller. Reproduced from [8]. To determine the net axial force, the pressure distribution is integrated over the respective surface areas on the suction and back side of the impeller. For a given pressure profile P (R), the axial force F acting on the impeller can be estimated using the following expression: F = 2π ∫ R1 Rs P (R) · R dR (2.8) where Rs is the radius of the pump shaft and R1 is the radius of the inner edge of the vane (as marked in Figure 2.3). From Figure 2.3, it is clear that the axial force that will act on the suction and the back side of the impeller are not equal, due to the difference in pressure distribution. To further illustrate, a schematic diagram of the axial forces acting on the semi-open impeller is shown in Figure 2.4. The net axial force T can be approximated by combining those individual contributions as: T = F1 − (F2 + F3 + F4) (2.9) 11 2. Theory where, F1 is the pressure force on the rear shroud (back side of the impeller), F2 is the opposing pressure force on the impeller hub, F3 is the force caused by the pressure differential across the blade tip clearance (the tiny gap between the blade tips and the stationary casing), and F4 is the dynamic reaction force resulting from the momentum change of the incoming fluid at the impeller inlet. This force imbalance leads to a net axial thrust directed along the shaft. In vertically oriented pumps, such unbalanced axial forces can exert a considerable amount of force on the bearing, potentially affecting the mechanical reliability of the pump. To mitigate this axial load, several balancing techniques, such as wear rings, back vanes, and balancing holes, are employed to redistribute pressure and reduce the load on the bearing. Notably, the axial force tends to decrease with increasing flow rate [10]. The following sections cover the underlying theory and function of balancing holes in a semi-open impeller. Pump Shaft Rear Shroud Impeller Vanes Figure 2.4: Schematic diagram of axial force on a semi-open impeller. 2.2.2 Balancing Holes Balancing holes (as illustrated in Figure 2.5) are typically small openings drilled through the back side (rear shroud) of the impeller in a centrifugal pump. These holes hydraulically connect the high-pressure region on the back side of the impeller to the low-pressure region near the suction side. By allowing a controlled reverse flow from the back side to the suction side, the balancing holes help equalize the pressure across the impeller, thereby significantly reducing the axial load [11]. This makes them a simple and effective method for reducing axial force in centrifugal pumps [2]. 12 2. Theory Balancing Hole Pump Shaft Rear Shroud Impeller Vanes Figure 2.5: Schematic showing the balancing hole in a semi–open centrifugal impeller. 2.2.2.1 Effect of Balancing Holes on Pressure Distribution and Axial Force As discussed in the previous section, the introduction of balancing holes in the impeller results in a reduced pressure distribution compared to the case without balancing holes (refer to Figures 2.3 and 2.4). This modification produces a less steep radial pressure gradient across the back shroud. This reduction in backside pressure lowers the axial force acting in that region, thereby reducing the net axial force acting on the impeller. This demonstrates the effectiveness of balancing holes in mitigating axial thrust and improving the mechanical balance of semi-open impellers. 2.2.2.2 Design Configurations While adding balancing holes is beneficial in mitigating the axial force, it might lead to performance loss (head and efficiency) due to the secondary leakage flow within the pump. The severity of this trade-off primarily depends on the geometric config- uration of the balancing hole, and to a lesser extent, on the operating condition of the centrifugal pump. Therefore, the design of balancing holes must be optimally configured to mitigate as much axial force as possible while minimizing the adverse effects on the pump performance. The parameters that mainly define the design configurations of balancing hole in- clude hole diameter, number of holes, chamfer (the beveled edge at hole entrance), radial and circumferential position. Each of these design parameters has a significant effect on both axial force and pump performance. 2.2.3 Clearance Area In centrifugal pumps with semi-open impellers, the clearance area refers to the nar- row radial gap between the rear shroud of the impeller and the stationary casing (bottom of the motor). This region, shown in Figure 2.6, forms a leakage path through which high-pressure fluid from the discharge side can move toward the 13 2. Theory backside of the impeller. When fluid leaks through this gap, it enters the back side of the impeller, affecting the pressure distribution in that region. The pressure on the back side of the impeller is determined by the balance of leakage flow and the dynamics of the rotating impeller and stationary casing, and it is typically between the suction and discharge pressures. The clearance area (AClearance), which governs the leakage flow rate and contributes to axial thrust, is given by the expression in Equation 2.10 AClearance = DMean · s (2.10) where DMean is the mean diameter of the impeller and s is the radial clearance between the rear shroud of the impeller and the stationary casing (bottom of the motor). Pump Shaft Rear Shroud Impeller Vanes Back Vane Clearance Bottom of Motor Figure 2.6: Clearance area on a semi-open impeller. This area plays a critical role in estimating internal leakage and is a key consideration in the design of other pump components, such as balancing holes, to ensure proper pressure distribution across the impeller to minimize axial force. 2.3 Computational Fluid Dynamics CFD is a branch of fluid mechanics that uses powerful numerical techniques to analyze and predict fluid flow behavior in complex engineering problems. The foun- dation of CFD lies in the fundamental governing equations of a fluid motion, which for isothermal incompressible flow with constant viscosity are the following: • Continuity equation (Conservation of Mass): ∂vi ∂xi = 0 (2.11) • Navier–Stokes equation (Conservation of Momentum): 14 2. Theory ρ ∂vi ∂t︸ ︷︷ ︸ Transient term + ρ ∂vivj ∂xj︸ ︷︷ ︸ Convective term = − ∂p ∂xi︸ ︷︷ ︸ Pressure gradient + µ ∂ ∂xj ( ∂vi ∂xj )︸ ︷︷ ︸ Viscous (Diffusion) term (2.12) where, vi and vj denote the components of the velocity vector, p is the pressure, ρ is the fluid density, and µ is the dynamic viscosity. These equations are discretized and solved iteratively to predict and analyze the dis- tribution of velocity and pressure within the fluid domain. In case of laminar flow, these equations are straightforward to solve. However, for turbulent flows—which are common in centrifugal pumps—directly solving these equations become much more complex due to their chaotic, nonlinear and multiscale nature. To make the problem solvable, the Reynolds-Averaged Navier–Stokes (RANS) equa- tions are used. In this method, the flow variables such as velocity and pressure are split into mean and fluctuating component. Here, the mean component represents time-averaged or ensemble-average value of the flow variable and the fluctuating component represents the deviation of the instantaneous value from the mean. vi = v̄i + v′ i (2.13) p = p̄ + p′ (2.14) v̄ = 1 2T ∫ T −T v dt (2.15) where vi is the instantaneous velocity component, v̄i is the time-averaged (mean) velocity, v′ i is the fluctuating velocity component, p is the instantaneous pressure, p̄ is the mean pressure, p′ is the pressure fluctuation, and v̄ is the time-averaged scalar velocity. Substituting Equation 2.13 & 2.14 into Equation 2.12 and performing time averaging (along with some mathematical manipulation) gives the Reynolds Averaged Navier- Stokes equation: ∂vi ∂xi = 0 (2.16) ρ ∂vi ∂t + ρ ∂v̄iv̄j ∂xj = − ∂p̄ ∂xi + ∂ ∂xj ( µ ∂v̄i ∂xj − ρv′ iv ′ j ) (2.17) In the resulting RANS momentum equation, a new additional term (τij = −ρv′ iv ′ j) has appeared in the diffusion term. This term is known as the Reynolds Stress Tensor, which accounts for the additional momentum transfer caused by turbulent fluctuations. This term represents the effect of turbulence on the mean flow and introduces 6 additional unknown terms, on top of the four original unknowns (3 velocity components and pressure) in the system of equations. To solve this problem, 15 2. Theory the additional term that appeared in Equation 2.17 must be modeled. The process of modeling those terms is known as turbulence modeling, which will be explained in the upcoming section. 2.3.1 Turbulence Modeling Turbulence modeling is a mathematical approach used in CFD to predict turbulent flow behavior within a fluid system. Due to its inherently chaotic, unstable, and fluctuating nature, accurately modeling this is a complex task. Instead of resolv- ing all the turbulent scales directly (Direct Numerical Simulation method), which is computationally expensive, a turbulence model has been employed to approximately represent their influence on the mean flow. Hence, different turbulence models have been developed to model the unknown additional term (Reynolds stress) and close the system of equations. In this study, the cubic k − ε turbulence model is used to analyze the internal flow effects within the centrifugal pump, due to its improved accuracy in capturing the anisotropic effects which is common in turbomachinery. The upcoming section provides a brief overview of standard k − ε turbulence model then discussing about cubic k − ε model. 2.3.2 Standard k − ε Turbulence Model The standard k − ϵ model is the most commonly used turbulence model to provide a mathematical representation of the turbulent flow field. It is a two-equation model, meaning it uses two additional transport equations, namely, turbulent kinetic energy (K) and turbulent dissipation rate (ϵ) along with the RANS equation, to help solve and predict the turbulent flow properties. The Boussinesq hypothesis models/relates the additional term (−v′ iv ′ j - Reynolds Stress Tensor) in the RANS equation to the mean velocity gradient using eddy (νt) viscosity, through a linear-relationship to close the RANS equation. The mathematical expression of the hypothesis is expressed as: −u′ iu ′ j = 2νtSij − 2 3kδij (2.18) Sij = 1 2 ( ∂vi ∂xj + ∂vj ∂xi ) (2.19) aij = −u′ iu ′ j + 2 3kδij = 2νtSij (2.20) νt = Cµ k2 ε (2.21) where Cµ is a constant with a value of 0.09, aij is the anisotropic part of the Reynolds stress tensor, δij is the Kronecker delta and Sij is the mean strain rate tensor. By applying the Boussinesq hypothesis, the 10 unknown terms in the RANS equations are reduced to a solvable system with 2 additional unknowns: the turbulent kinetic 16 2. Theory energy and the turbulent dissipation rate. These can be determined by solving their respective transport equations. However, this linear model as said earlier assumes an isotropic (which means the turbulence has the same intensity and statistical prop- erties in all 3 directions) relationship between the turbulent stresses and the mean strain rate, which limits its ability to accurately predict complex flows with strong anisotropy, streamline curvature, or separation. The transport equation for the turbulent kinetic energy: ∂(ρk) ∂t︸ ︷︷ ︸ Unsteady term + ∂(ρkvj) ∂xj︸ ︷︷ ︸ Convection = ∂ ∂xj [( µ + µt σk ) ∂k ∂xj ] ︸ ︷︷ ︸ Diffusion + Pk︸︷︷︸ Production − ρε︸︷︷︸ Dissipation (2.22) The transport equation for the turbulent dissipation rate: ∂(ρε) ∂t︸ ︷︷ ︸ Unsteady term + ∂(ρεvj) ∂xj︸ ︷︷ ︸ Convection = ∂ ∂xj [( µ + µt σε ) ∂ε ∂xj ] ︸ ︷︷ ︸ Diffusion + C1ε ε k Pk︸ ︷︷ ︸ Production − C2ερ ε2 k︸ ︷︷ ︸ Destruction (2.23) Where C1ε, C2ε, and C3 are model coefficients whose values depend on the specific form of the k-ε model chosen. 2.3.3 Cubic k − ε Turbulence Model In this model, the Reynolds stress is modeled using a more general non-linear relation that includes both the mean strain rate tensor and the mean rotation rate tensor, also known as the vorticity tensor [12]. This general form is expressed as: aij = ∑ n GnT (n) ij (2.24) aij = − 2Cµτ Ŝij + c1τ 2 ( ŜikŜkj − 1 3 ŜℓkŜkℓδij ) + c2τ 2 ( Ω̂ikŜkj − ŜikΩ̂kj ) + c3τ 2 ( Ω̂ikΩ̂jk − 1 3Ω̂ℓkΩ̂kℓδij ) + c4τ 3 ( ŜikŜkℓΩ̂ℓj − Ω̂iℓŜℓkŜkj ) + c5τ 3 ( Ω̂iℓΩ̂ℓmŜmj + ŜiℓΩ̂ℓmΩ̂mj − 2 3Ω̂mnΩ̂nℓŜℓmδij ) + c6τ 3ŜkℓŜkℓŜij + c7τ 3Ω̂kℓΩ̂kℓŜij (2.25) Ω̂ij = 1 2 ( ∂vi ∂xj − ∂vj ∂xi ) (2.26) where T (n) ij are the basis tensors and Gn are scalar coefficients. This nonlinear Boussi- nesq approximation allows the model to more accurately capture anisotropic effects, secondary flows, and curvature-related phenomena, which are especially important for simulating the complex turbulent flows in centrifugal pumps. 17 2. Theory 2.3.4 Multiple Reference Frame model - Frozen Rotor Ap- proach The Multiple Reference Frame model is the most widely used CFD approach for simulating rotating machinery like turbines, fans and pumps. It is a steady-state simulation method where the entire computational domain is divided into multiple zones – a rotating zone (typically an impeller or rotor) and a stationary zone (volute or casing). The rotating and stationary zones are mathematically coupled through a non-sliding interface (fixed interface) where the flow variables are interpolated be- tween the frames of reference. Frozen Rotor Approach is a specific implementation of the MRF model, where the rotating region is held fixed at a certain angular position with respect to the sta- tionary region. The governing equations are solved in a rotating reference frame for the rotating region and a stationary reference frame for the stationary region. Even though the rotating region does not physically rotate during the simulation, the rotational effects are fully accounted for in the governing equations by additional source terms, especially centrifugal and Coriolis forces in the momentum equations. Stationary Zone Moving Reference Frame Zone Interface Figure 2.7: Geometry with one rotating impeller. Taking an example of a mixing tank with a single impeller as illustrated in Figure 2.7, the impeller along with the surrounding flow (orange region) is defined to be the moving frame, and the flow outside the impeller is defined as the stationary frame. The interface between the two frames is assumed to have a steady-state condition, where it serves as a surface of revolution for the two frames [13]. The interface (non-sliding) between the two zones allows momentum and pressure to be transferred without the need for any transient time-stepping scheme. The Frozen Rotor Approach provides a computationally effective way to evaluate the performance of rotating machines. It is suitable for steady-state performance pre- dictions, initial design and performance assessment and a precursor to more detailed transient CFD analysis. However, this approach is not applicable for simulations in- volving transient effects like unsteady interactions between the rotor and the stator, 18 2. Theory blade passing phenomenon, and flow-induced fluctuations that require time-accurate modeling. 2.4 Multi-Objective Optimization Multi-objective optimization refers to an optimization problem that involves more than a single objective function, which needs to be optimized simultaneously, often with conflicting goals. In such problems, the aim is to find solutions that represent the best possible trade-offs between those objectives, since one objective may lead to the degradation of the other objective. These solutions are simply called Pareto- optimal solutions or non-dominated solutions. 2.4.1 Pareto Optimality In multi-objective optimization problems, it is possible to obtain multiple feasible solutions rather than a single optimal solution. Identifying a single feasible solution in such optimization processes is quite challenging due to the inherently conflicting nature of each objective. Instead, a more meaningful strategy is to determine a set of feasible solutions that represent the best possible trade-offs, a concept known as Pareto optimality. A solution is said to be Pareto optimal if no other solution in the objective space ex- ists that improves one objective without worsening the other objective. These types of solutions available in the design / objective space are also called non-dominated solutions. These non-dominated solutions collectively form the Pareto front, which defines the boundary of optimal trade-offs in the objective space. This is defined more precisely as: A point x∗ in the feasible design space S is Pareto optimal if and only if there does not exist another point x in the set S such that f(x) ≤ f(x∗) with at least one fi(x) < fi(x∗) [14]. The collection of such non-dominated solutions can be visualized as forming a Pareto front in the objective space. In problems with two objectives, as shown in Figure 2.8, the Pareto front typically appears as a smooth curve connecting the non-dominated solutions. Thus, the Pareto optimal set is always on the boundary of the feasible objective space. Each point on the Pareto front represents a valid trade-off solu- tion, and no point is strictly better than another across all objectives. Therefore, the final selection among Pareto-optimal solutions depends on the decision maker’s preferences and contextual priorities. While the Pareto-optimal set is always on the boundary of the feasible objective space in a constrained scenario, it is significant to note that this boundary is not always defined by the explicit constraints. The Pareto front also emerges in un- constrained problems, where the optimal trade-offs purely arise from the inherent 19 2. Theory conflict between objectives. In such problems, the Pareto optimal set is determined by comparing the non-dominated solutions across the objective space or analyti- cally by identifying the points where the gradients of the objective functions oppose each other, indicating that no common descent direction exists for simultaneous improvement. Pareto Optimal Points Figure 2.8: Pareto front with non-dominated optimal points. 2.4.1.1 Preferences and Utility Function Mathematically, there are infinitely many Pareto-optimal solutions available within an objective space. This makes the multi-objective optimization problem more chal- lenging to solve and find an optimal solution. Apparently, it becomes necessary to make decisions regarding which solution is preferred. This can be done by the process of providing preference to the objective function. Preferences reflect the decision maker’s priorities or opinions concerning different solutions within the ob- jective space. Ideally, a multi-objective method should incorporate these preferences so that the outcome better aligns with the decision maker’s needs. However, rep- resenting preferences mathematically is a tedious and non-trivial task, especially when aiming for a multi-objective target function that covers all the nuances of the decision maker’s need. However, various methods have been developed to incorporate the preferences math- ematically into an objective function. Among those, the most traditional and math- ematically straightforward approach is to assign weights to each objective based on their relative importance, allowing them to reflect the decision maker’s preference [14]. Although straightforward, it is not an optimal strategy, as it may cause the optimization process to overly focus on a single objective, completely deviating from the core principle of a multi-objective optimization. To tackle this, Steuer [15] pro- vided a mathematical formulation using a hyperplane for interpreting and assigning weights in a multi-objective target function. Based on this, the mathematical for- 20 2. Theory mulation for calculating weights in a multi-objective optimization process can be expressed as: ∆i = 1 N N∑ j=1 |fi(xj) − fi,ref| (2.27) where fi(xj) is the value of the ith objective for the jth sample, fi,ref is the reference value for the ith objective, and N is the total number of samples. Equation 2.27 represents the mean absolute deviation from its reference value. Then, the inverse of the mean absolute deviation is computed to emphasize objectives that are closer to their reference values: Ii = 1 ∆i (2.28) The normalized weights wi are obtained by: wi = Ii∑m k=1 Ik (2.29) where wi is the normalized weight assigned to the ith objective, Ii represents the inverse of the mean absolute deviation for the ith objective, Ik denotes the inverse deviation of the kth objective (with k ranging over all objectives), and m is the total number of objectives considered. While the weightage method represents the preferences among different objectives, the utility function plays the role of combining these multi-objective functions into a single-objective framework (i.e., by scalarizing multiple objectives into a single objective function), thereby simplifying the optimization process. Scalarization is the mathematical process of combining or transforming a set of different objective optimization problems into a single objective function [14]. By scalarizing, the op- timization process becomes easier, and it provides a structured and accurate way to analyze the trade-offs between different sets of objectives in multi-objective op- timization. The general form of the scalarized objective function is expressed as shown in Equation 2.30. F (x) = m∑ i=1 fi(x) (2.30) where F (x) is the scalarized objective function, fi(x) is the ith objective function, and m is the total number of objectives. This transformation allows the selection and evaluation of the optimal solution to be carried out more mathematically without neglecting the decision maker’s demands while properly taking into account the trade-offs between different objectives within a unified approach. 2.4.2 Normalization Normalization is the process of scaling or transforming a set of absolute data sets relative to a reference value, typically to make the data sets into dimensionless 21 2. Theory or comparable form, allowing consistent comparison across magnitudes, units, or ranges. Normalization is widely used to omit the effects of scale, making analyses based on relative differences rather than absolute values. In this study, different normalization technique were used to quantify the deviation between simulated results and a reference data. Below is a brief explanation of the techniques used in this study. 2.4.2.1 Relative Squared Error (RSE) The Relative Squared Error (RSE) is a widely used metric system to analyze or quantify how close the observed (in this case, simulated result) value deviates or matches with the reference value. It is simply the ratio of the total deviation of all observed values from the reference values to the total deviation of the same observed values from their mean values. By this, the RSE brings the deviation to a relative scale. The general form of RSE, often used for assessing the performance across an entire dataset, is expressed as: RSE = ∑n i=1(yref,i − ŷi)2∑n i=1(yref,i − yref )2 (2.31) where yref,i is the reference value for the ith data point, ŷi is the observed or simu- lated value for the ith data point, yref is the mean of all reference values, and n is the total number of data points. In this study, a modified version of RSE, known as the Pointwise Relative Squared Error, has been used to normalize and evaluate the deviation of individual data points with respect to the reference values. 2.4.2.2 Pointwise Relative Squared Error The modified version (for this study) known as the Pointwise RSE focusses on ana- lyzing the deviation of each individual dataset from the reference value. This offers a more detailed understanding of the performance at each specific observation, rather than presenting a single overall measure of error of the observed data. The Pointwise RSE for an individual dataset is expressed as: Pointwise RSE = (yref,i − yi)2 (yref, i)2 (2.32) where yi is the observed or simulated value for the ith configuration, and yref, i is the corresponding reference value for the ith configuration. The necessity and the practical implementation of this technique are discussed in Section 3.4.1. 22 2. Theory 2.4.3 Bayesian Optimization Bayesian optimization is a technique used for finding optimal solutions for objective functions that are expensive to evaluate since they provide no information about their gradients [16]. This function can be noisy or non-convex in nature further complicating the traditional optimization like the first order and second order based methods. Instead of directly evaluating and optimizing the true objective function, which can be computationally expensive and time-consuming, Bayesian optimiza- tion relies on constructing a surrogate model to approximate the function and guide the search for the optimal solution of the objective function. In Bayesian Optimization, the Gaussian process as the surrogate model provides a probabilistic prediction of the true black-box objective function value and also an estimate of the uncertainty at any point in the design space. This dual information associated with an acquisition function enables mathematical decision-making on where to explore or sample next based on the current information by balancing ex- ploration (where the model is uncertain) and exploitation (where the model predicts high performance) [16]. 2.4.3.1 A Bayesian Approach to Optimization The Bayesian optimization is fundamentally concerned with finding the minimizer of the objective function without needing the access to the Gradient or Hessian of the function [16]. Rather, it queries on the function f at a given point x and re- turning the value of f(x). In the Bayesian approach, instead of attempting to find the optimal point x for a function f , it aims to learn a best model of the function based on previously observed points x. By conditioning on these observed points, the Bayesian optimization extracts an expected function f and a confidence interval, which are used to guide the selection of future sampling points. This makes Bayesian optimization particularly well suited for complex design opti- mization and computationally expensive simulation tasks. In this study, the method is employed to optimize different balancing hole configurations in a centrifugal pump, where each design evaluation involves computationally intensive CFD simulations. The following sections describe the two core components of Bayesian optimization: the Gaussian Process and the acquisition function. 2.4.3.2 Gaussian Process A Gaussian process is a probabilistic non-parametric model used in Bayesian opti- mization to construct a surrogate model over the objective function. It defines a probability distribution over the functions and can be updated with new observations to refine its predictions [16]. The Gaussian process provides not only an estimate of the objective function’s value at any input points but also an uncertainty measure at those same points in the input space. 23 2. Theory f(x) ∼ GP(µ(x), K(x, y)) (2.33) where µ(x) is the mean function, assigning to each input point x the expected value E[f(x)], and K(x, y) is the covariance function (or kernel), assigning to each pair of input points x and y the covariance between f(x) and f(y). The covariance function is also called the kernel of the Gaussian process. It simply measures the similarity or correlation between the values of f(x) and f(x′). Among various kernel functions available, the most widely used kernel in Gaussian process is the Radial Basis Function (RBF) [16] kernel also known as the squared exponential kernel or Gaussian kernel, defined as: K(x, y) := k · exp { −∥x − y∥2 2 2σ2 } , (2.34) where k is the variance parameter that sets the overall scale of variation, and σ is the length-scale parameter that controls the smoothness of the function. This kernel function ensures that the input closer in the design space have highly correlated function values, which distant inputs have low level of correlation. By using this kernel function, Gaussian process surrogate model captures the structure of the underlying true objective function effectively, allowing Bayesian Optimization to propose new promising candidates [17]. 2.4.3.3 Acquisition Function The acquisition function is a guideline that suggests where to query the next sample points based on the posterior predictive distribution (refer to Section 2.4.5) of the function provided by the Gaussian process. One of the most commonly used ac- quisition functions is the expected improvement (EI), which quantifies the expected amount by which a new sample will improve the best observed objective value [16]. Let f ∗ be the best value observed so far, and let f̃(x) ∼ N (µ(x), σ2(x)) be the posterior predictive distribution of the Gaussian Process at input x. The expected improvement (EI) is defined as: EI(x) = E[max(0, f ∗ − f̃(x))] = (f ∗ − µ(x)) · Φ(Z) + σ(x) · ϕ(Z) (2.35) Z = f ∗ − µ(x) σ(x) (2.36) where Z is the standardized improvement, Φ(Z) is the cumulative distribution func- tion (CDF) of the standard normal distribution, and ϕ(Z) is the probability density function (PDF) of the standard normal distribution. The next sampling point is chosen by maximizing the acquisition function defined in Equation 2.35, where the goal is to identify the location in the design space that 24 2. Theory offers the greatest expected improvement over the current best observation. This formulation balances exploitation and exploration, making expected improvement (EI) a powerful and practical choice for design optimization tasks. 2.4.3.4 Bayesian Optimization Work Flow Start Initial Sampling (Sobol sequence sampling) Fit Gaussian Process surrogate model Maximize Acquisition Function (EI) Evaluate Objective Function at Selected Point Update Model with New Data Convergence Met? End Yes No Figure 2.9: Bayesian optimization workflow. Figure 2.9 illustrates the overall process involved in the Bayesian optimization work- flow. It can be seen from the flow chart that the workflow starts with an initial sampling step using Sobol sequence. This initial step ensures that the optimization starts with a set of input points that are uniformly (evenly) distributed in the design space, providing a solid foundation for Bayesian Optimization to build the surrogate model using Gaussian process. The theoretical background for the Sobol sequence sampling is presented in the following section. 25 2. Theory 2.4.4 Sobol Sequence Sampling In high-dimensional design and optimization problems, especially those involving surrogate modeling or Bayesian optimization, the initial sampling of design points plays a critical role in determining the efficiency and accuracy of the optimization process. Traditional random sampling often leads to clustering and poor coverage of the design space. To overcome this, low-discrepancy sequences, such as the Sobol sequence, are employed to provide better uniformity and space-filling properties. A Sobol sequence is a type of quasi-random sequence designed to fill the unit hyper- cube [0, 1]d as uniformly as possible. Unlike purely random (Monte Carlo) samples, which suffer from uneven spacing and gaps in the domain, Sobol sequences exhibit low discrepancy, meaning the generated points are distributed uniformly within the design space [18]. This makes them highly effective in areas such as numerical in- tegration, global sensitivity analysis, and initial design of experiments (DoE) for machine learning and optimization. Sobol’ sequences are a type of digital (t, s)- sequence in base 2, according to the definition or framework given by Niederreiter [19]. 2.4.4.1 Discrepancy and Uniformity Properties Sobol sequences are constructed to minimize discrepancy—a measure of the devi- ation of a sequence from a perfectly uniform distribution. The discrepancy DN of the first N Sobol points satisfies: DN = O ( (log N)d N ) This is asymptotically superior to the convergence rate of Monte Carlo sampling: DN = O(N−1/2) 2.4.5 Posterior Predictive Distribution The posterior predictive distribution, a fundamental concept in Bayesian statistics, is the distribution of possible unobserved outcomes based on the Bayesian model [20]. It is obtained by integrating the entire posterior distribution of the model parameters by averaging over the likelihood of all possible parameter values weighted by their posterior probability. The posterior predictive distribution for a new input x∗ is given by: p(y∗ | x∗, D) = ∫ p(y∗ | x∗, θ) p(θ | D) dθ (2.37) where θ denotes the model parameters, D is the observed dataset, and y∗ is the predicted output at the new input x∗. 26 3 Methodology 3.1 Design configurations for balancing hole opti- mization To develop a robust multi-objective target function for the design optimization of balancing hole, a set of well-defined suitable design parameters were considered in this study. This set of design parameters includes the diameter of the balancing hole, the number of holes, their radial position relative to the center of the impeller, and their circumferential position relative to the impeller’s front and rear vanes. The selection of these design parameters was highly influenced by their impact on axial force reduction, leakage flow management and overall hydraulic performance of the pump. (a) Pressure / Back side. (b) Suction / Front side. Figure 3.1: Semi-open impeller without balancing hole. The 3D CAD model of the impeller along with the balancing hole layout, was pro- vided by Xylem. The geometrical representation of the working impeller without balancing hole is shown in Figure 3.1. The design configurations explored in this study were defined using the Relation feature (an inbuilt option) within PTC Creo, allowing for a systematic control over the diameter, radial and circumferential po- sitioning of the balancing holes. These parametric relations were predefined by Xylem and reused in this study to design multiple balancing hole configurations. 27 3. Methodology The following subsections describe how each of these parameters was defined. 3.1.1 Number of balancing holes The number of balancing holes used in this optimization study was fixed at 3. This is because the impeller has three blade passages, and each hole is placed to match one passage. So, there is one hole for each blade passage, as shown in Figure 3.3a. 3.1.2 Balancing hole diameter The diameter of the balancing hole range was determined using the clearance area (AClearance), which was provided by Xylem. According to the guideline proposed in Gülich’s centrifugal pump design framework [3], the total area of the balancing hole should be 4–5 times the clearance area. Figure 3.2, illustrates red arrow highlights the clearance between the impeller and the station casing (bottom of motor), while the blue arrow indicates the washout channel designed to remove debris entering through the clearance. Together, these two regions constitute the total clearance area used for balancing hole sizing. Figure 3.2: Clearance area of the semi-open impeller. Given three balancing holes fixed, the corresponding minimum and maximum range of total area were calculated as: ATotal ∈ [1 × AClearance, 5 × AClearance] Aper Hole = [1 × AClearance 3 , 5 × AClearance 3 ] = [Aper Clearance,min, Aper Clearance,max] Given that the holes are circular, the diameter range is computed as: 28 3. Methodology D = √ 4 × Aper Hole π ⇒ √4 × Aper Hole,min π , √ 4 × Aper Hole,max π  = [Dmin, Dmax] Thus, the balancing hole diameter considered in this study ranges from Dmin and Dmax. 3.1.3 Circumferential position of balancing hole Defining the circumferential positioning of the holes required a two-part approach due to the geometry of the impeller. The impeller was conceptually divided into the back shroud, and the front blade section. Drilling through the impeller based solely on the back-side layout would risk interference with the front blades. To address this, the front section of the impeller was rotated about its central axis to align the front and rear hole positions as illustrated in Figure 3.3b. (a) Number of holes. (b) Rotation in front side. Figure 3.3: Semi-open impeller with balancing holes. 3.1.3.1 Back Side of the Impeller To define hole positions on the rear shroud, a geometric method based on back vane spacing was adopted. A centerline was drawn along one of the back vanes and ro- tated around the impeller’s central axis. Since the back vanes are uniformly spaced, the angular separation between two vanes is 60◦. Therefore, the circumferential positioning of holes will be explored between the two adjacent vanes as shown in Figure 3.4a 3.1.4 Radial Position of Balancing Hole The radial position of the balancing hole was defined using two diameters, D1 and D2, as shown in Figure 3.4a. These diameters represent the radial limit within which the hole positioning was explored. Within this bounded radial region, a curved guide 29 3. Methodology path—shown as the red arc in Figure 3.4a—was drawn to follow the same curvature used for the circumferential placement of the holes. This arc is similar to the shape of the back vane and maximizes the available surface area for hole placement. This curvature-based method ensures that the hole positioning effectively explores the available space between the adjacent vanes (the red zone in Figure 3.4b), enabling a broader coverage of the impeller rear shroud region. The radial distance from the impeller center to each hole was then systematically varied along this curved path during the design exploration process, within a range of 0.203 to 0.342. Here these two values represent the normalized value of D1 and D2 with respect to the impeller diameter. (a) Radial position. (b) Feasible region for hole placement. Figure 3.4: Balancing hole placement and feasible region in the semi-open impeller. 3.1.5 Finalization of Design Parameter Ranges Although the diameter, radial, and circumferential position of the balancing holes were initially defined independently, they were ultimately coupled during the opti- mization process. In practice, this coupling introduced new geometrical constraints that necessitated refinement of the original parameter ranges. Specifically, the cir- cumferential angle range between 0◦ and 60◦ was excluded, as the placements within this range would intersect with the back vanes, potentially compromising the vane structure. Similarly, the radial positioning range was adjusted to prevent the hole position from being too close to either the impeller center or its peripheral edge. In addition to this, the varying range of hole diameter contributed to these constraints, as larger diameters increased the risk of interference with key structural features. Furthermore, the front-side rotation required for through-hole drilling was adjusted to ensure geometric feasibility across all configurations. To effectively explore the geometric space while avoiding interference with critical impeller features, the cir- cumferential and radial position ranges were carefully refined to avoid interference with any structural features within the impeller geometry. The final constrained ranges used in this study for the design parameters are summarized in Table 3.1. 30 3. Methodology Table 3.1: Design parameter ranges before and after geometric refinement. Design Parameter Initial Range Constrained Range Circumferential Position 60◦ 24◦ Radial Position normalized with Impeller Diameter 0.203 to 0.342 0.245 to 0.3 Rotation 120◦ 66◦ 3.2 Computational Modeling This section outlines the computational modeling strategy used for simulating the internal flow within the centrifugal pump considered in this study. Both the sim- ulation setup and meshing strategies were provided by Xylem. While a base- line mesh for the entire computational domain was provided , local mesh refinement was performed near the balancing hole regions to re- solve the complex flow phenom- ena more accurately. All simulations were conducted using CFD++, with mesh generation and pre-processing performed in ANSA. 3.2.1 Computational Domain The computational domain, as shown in Figure 3.5, represents the internal geometry of the centrifugal pump used in this study. The Frozen Rotor method was employed to model the rotor and stator interactions within the computational domain. This approach minimizes computational cost and time compared to transient simulations. To facilitate this modeling strategy, the computational domain is divided into two regions (refer to Figure 3.5): Rotating Region Stationary Region Figure 3.5: Computational domain. 31 3. Methodology • Rotating Region: The rotating region represents the part of the computa- tional domain that is subjected to rotational motion. It was assigned a fixed rotation speed of 990 RPM to replicate actual pump operating conditions. Flow within this region is solved in a rotating frame of reference. However, it can be seen from Figure 3.5 that certain boundaries included within this region do not physically rotate in reality. This will be summarized in the up- coming Section 3.2.2. The computational domain is extended upstream, as shown in Figure 3.5, to include a stationary inflow section that minimizes the influence of numerical boundary conditions on the physical flow development near the impeller. Additionally, the rotating region is extended slightly up- stream of the impeller to accurately capture the incoming flow characteristics and to provide a sufficient buffer for the smooth interpolation of flow variables across the rotating–stationary interface. This transition zone, necessitated by the use of different reference frames, helps maintain continuity in flow variable gradients, reduces numerical artifacts, and ensures accurate flow prediction at the impeller inlet. • Stationary Region: This region represents the stationary part of the compu- tational domain. The components assigned to this region remain fixed during the simulation and interact with the rotating region through interface bound- ary. Although no volute is included in the current model and the domain appears ax- isymmetric, the stationary region is retained to allow for potential inclusion of a volute in future studies, which would break axisymmetry and require a fixed region. 3.2.2 Boundary Conditions To ensure accurate flow modeling, appropriate boundary conditions were assigned to each surface of the computational domain. The conditions were assigned based on the physical behaviour of each component and the simulation strategy. Since the Frozen Rotor simulation strategy were used to handle the interaction between the rotating and stationary regions, boundaries in those regions were defined either in the global (stationary) frame of reference or local (rotating) frame of reference, depending on the actual motion of the components. As shown in Figure 3.6, the boundaries were labeled according to the Property IDs (PIDs) in the ANSA pre-processing environment. The names are consistent with the physical surfaces represented in the geometry and distinguish between rotating and stationary regions. The assigned condition for each boundary is summarized as: 1. Outlet: Pressure Outlet. 2. Rotating Smooth: Wall (Viscous, Wall Function). 3. Slip: Wall (Inviscid). 32 3. Methodology 4. Inflow: Mass Flow Rate (Value depends on Operating condition of the pump). 5. Stat Rough In Rot: Wall (Viscous, Wall Function, Rotation about Z–axis, Rotation Rate: 0, Roughness: 0.00011 m). 6. Rotating Rough Back: Wall (Viscous, Wall Function, Stationary wrt Mesh Motion, Roughness: 0.00011 m). 7. Rotating Rough Channel: Wall (Viscous, Wall Function, Stationary wrt Mesh Motion, Roughness: 0.00005 m). 8. Stationary Smooth In Rot: Wall (Viscous, Wall Function, Rotation about Z–axis, Rotation Rate: 0). 9. Diffuser Slip: Wall (Inviscid). 10. Inlet Interface: Zonal (Patched Only, Simple Flow Through). 1 2, 6 and 7 10 3 8 9 5 4 Figure 3.6: Boundary condition for the computational domain. 33 3. Methodology Some boundaries located within the rotating region are assigned zero rotation to reflect the real case where only the impeller rotates. The surface roughness values were provided by Xylem and are based on experimental characterization of the pump component’s surfaces. The boundary condition types and recommended settings are based on the CFD++ user manual [21]. 3.2.3 Computational Mesh As previously mentioned in Section 3.2, the generated computational mesh in ANSA was provided by Xylem and validated through transient CFD simulations and mea- sured performance. Since the baseline mesh did not include balancing holes, lo- cal refinement was applied in those regions. This section first details the baseline meshing strategy for the original geometry, followed by a mesh independence study focused on the balancing hole. 3.2.3.1 Baseline Mesh Configuration (Without Balancing Hole) Figure 3.7 shows the sectional view of the computational domain with surface and volume meshes. The surface mesh uses a combination of quadrilateral and triangu- lar elements to effectively mesh complex areas like the impeller blades, back vanes, narrow passages, and sharp trailing edges. This hybrid approach ensures a good mesh quality with reduced skewness in complex regions and enables the creation of high-quality surface discretization. For the volume mesh, a combination of hexahe- dral, tetrahedral, and pyramidal elements was used throughout the computational domain. These elements are chosen for their favorable numerical properties, where hexahedral elements offer improved accuracy, reduced numerical diffusion, and en- hanced solver stability due to their structured nature, while tetrahedral and pyra- midal elements provide flexibility for meshing complex geometries. This meshing strategy helps in better resolution of the pressure and velocity gradients in intricate regions of a centrifugal pump geometry. Additional mesh details near the impeller and boundary layer regions are shown in Figures 3.8 and 3.9. The mesh settings are summarized in Table 3.2: Table 3.2: Mesh settings used in the original configuration Parameter Description / Value Surface Mesh Elements Hybrid mesh (Quadrilateral + Triangular) Surface Element Size Ranges from 0.8 mm to 30 mm Surface Mesh Growth Rate 1.2 Volume Mesh Type Hybrid mesh (Hexa + Tetra + Pyramid) Volume Element Size Maximum element size of the surface mesh Volume Mesh Growth Rate 1.2 Number of Boundary Layers 7 First Cell Height 0.2 mm Boundary Layer Growth Rate 1.2 Total Element Count 13.5 million cells 34 3. Methodology (a) Surface mesh. (b) Volume mesh. Figure 3.7: Computational mesh of the domain Figure 3.8: Boundary layer mesh. (a) Pressure / Back side. (b) Suction / Front side. Figure 3.9: Surface mesh on the semi-open impeller without balancing hole. 35 3. Methodology (a) Coarse mesh. (b) Fine mesh. Figure 3.10: Surface mesh refinement – back side (Balancing hole diameter = Dmax). (a) Coarse mesh. (b) Fine mesh. Figure 3.11: Surface mesh refinement – inner side (Balancing hole diameter = Dmax). 3.2.3.2 Mesh Refinement Study around Balancing Hole To ensure accurate resolution of flow characteristics near the balancing hole, a mesh independence study was carried out with a focus on localized refinement. Initially, surface mesh refinement was applied around the balancing hole on the impeller while maintaining the baseline volume mesh settings. Element sizes ranged from 0.8 mm to 7 mm. Figures 3.10 and 3.11 presents surface mesh examples for coarse and fine refinements, applied to the maximum diameter (Dmax) range of the balancing hole. 36 3. Methodology (a) Cylindrical shape. (b) Conical shape. Figure 3.12: Volume mesh refinement (Balancing hole diameter = Dmin). Due to the expected strong pressure gradients and jet-like flow behavior caused by throttling at the hole, volumetric mesh refinement was also tested. Two refinement zone geometries were evaluated.a cylindrical shape and a conical shape (refer to Fig- ure 3.12). These refinements increased local mesh density and improved resolution of the flow structures within and around the balancing hole. The final mesh settings selected through this study provided sufficient resolution for all tested balancing hole configurations and were used for simulation results discussed in Section 4.2. This mesh refinement study was performed at the BEP and included surface refinement for both the smallest (Dmin) and largest (Dmax) balancing hole diameters, while volume refinement was evaluated only for the smallest diameter configuration using a surface mesh element size of 1 mm. 3.2.4 Turbulence Modeling and Solver Settings The CFD simulations in this study were performed using a steady-state approach with solver settings and turbulence model configurations based on Xylem’s indus- trial CFD++ setup. A non-linear cubic k-epsilon model (refer to Section 2.3.3) was 37 3. Methodology used to model the turbulence. This model was chosen to accurately resolve the com- plex flow features that arise in centrifugal pumps with balancing holes, particularly the secondary flow caused by rotational effects and the related turbulent anisotropy, flow separation, and recirculation near both the balancing hole and the impeller blades. In addition to this, a non-equilibrium Launder-Spalding wall function was used to ensure accurate near-wall modeling in regions affected by strong pressure gradients and boundary layer separation, which are dominant near impeller blades and balancing holes. To maintain physical consistency in the turbulence quantities throughout the domain, the Schwartz realizability criterion was also selected. The numerical solver used an implicit time-stepping scheme with local time- step control to accelerate convergence while maintaining numerical stability. A second- order upwind scheme was used for spatial discretization to improve numerical ac- curacy and effectively resolve the flow gradients. Further, under-relaxation factors were used for both the coupled pressure-velocity solver and turbulence equations to ensure stability during iterations. The solver and the turbulence model settings used in this study are summarized in Table 3.3. These settings ensured stable conver- gence, accurate turbulence prediction, and robust performance across all operating conditions. Table 3.3: Turbulence modeling and solver settings Setting Value / Description Solver Type Pseudo-Steady-State Time Integration Local time stepping (Implicit) Turbulence Model Cubic k–ε Wall Function Methodology Launder–Spalding Wall Function Type Non-Equilibrium Wall Function Pressure Gradient Term Enabled Discretization Second Order Under-Relaxation Factor 0.5 (Coupled Solver), 0.95 (Turbulence) Working Fluid Incompressible water at 293.15 K 3.3 Initial Sampling for the Optimization Process using Sobol Sequence To initiate the optimization process and to ensure comprehensive coverage of the design space, the Sobol sequence was employed for sampling. As previoiusly men- tioned, four design parameters were considered in this study: balancing hole diame- ter (BH_Diameter), the circumferential position (BH_Circumferential_Back & BH_Circumferential_Front) in the front side of the impeller and back side of the impeller as well as its radial position (BH_Radial_Position). 38 3. Methodology 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 BH_Circumferential_Back B H _ D ia m et er 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 BH_Circumferential_Back B H _ C ir cu m fe re nt ia l_ Fr on t 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 BH_Circumferential_Back B H _ R ad ia l_ P os it io n 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 BH_Diameter B H _ C ir cu m fe re nt ia l_ Fr on t 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 BH_Diameter B H _ R ad ia l_ P os it io n 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 BH_Circumferential_Front B H _ R ad ia l_ P os it io n Figure 3.13: Pairwise Scatter Plots of Initial Sobol Samples. Due to the high computational cost associated with CFD simulations, 12 initial de- sign points were considered and it corresponds to 3 times the total number of design parameters. This sampling strategy aimed to strike a balance between adequate ex- ploration of the design space and computational efficiency, particularly in the early phases of the optimization. Furthermore, it ensured that the surrogate model was initialized with sufficient data to support the subsequent Bayesian optimization pro- cess. The 12 initial seed samples, each representing a unique design configuration, were generated using the Python code script provided in appendix A.1. The actual code for the Sobol sequence samples, developed by Xylem, is not included in this re- 39 3. Methodology port due to confidentiality restrictions. Figure 3.13 illustrates how the initial seed samples are distributed across the design space. Table 3.4 summarizes the nor- malized design parameter values—such as back, diameter, front, and radial loca- tions—for a subset of these Sobol-generated configurations. The performance re- sults of these 12 design configurations, evaluated at four different flow rates (Q∗ = 0.50, 0.75, 1.00, and 1.25), are discussed in detail in the results chapter under section 4.8. Table 3.4: Design parameters and their spatial positions for the 12 initial Sobol-seeded configurations. Configuration Back position Hole diameter Front position Radial position 1 0.678 0.858 0.903 0.719 2 0.038 0.067 0.166 0.007 3 0.333 0.534 0.705 0.885 4 0.944 0.258 0.480 0.330 5 0.781 0.704 0.118 0.443 6 0.426 0.495 0.826 0.778 7 0.188 0.904 0.304 0.142 8 0.581 0.179 0.504 0.571 9 0.512 0.618 0.329 0.872 10 0.156 0.331 0.598 0.411 11 0.481 0.802 0.031 0.541 12 0.874 0.030 0.792 0.236 3.4 Target Function Development Strategy Figure 3.14: Workflow for multi-objective target function development. To construct a robust multi-objective target function for the optimization of balanc- ing hole design configurations, a systematic approach was employed. This approach involved the integration of key performance metrics, such as head developed, power consumption, and axial force, into a unified evaluation framework. The resulting target function enables the rankings of design configurations based on their ability to maintain hydraulic performance while minimizing axial force. The development and 40 3. Methodology refinement process of the target function is outlined in Figure 3.14, which highlights the key steps from the evaluation of the initial samples generated by the Sobol se- quence to function refinement. This section describes the approach under two main categories: normalization strategy and target function formulation. 3.4.1 Normalization Strategy Normalization procedure is a crucial component of this study, since the considered performance parameters have different units and magnitudes. Without appropriate normalization, the target function could be biased toward parameters with larger numerical values, potentially undermining the optimization process. The normal- ization procedures adopted are discussed in the following subsections. 3.4.1.1 Normalization of Head and Power Initial normalization was carried out using the standard deviation derived from the initial design configurations data obtained using the Sobol sequence sample. Al- though statistically valid, this approach introduced challenges due to the uneven spread of normalized values, which made it challenging to construct a stable and in- terpretable target function. The normalization formula using the standard deviation is expressed as: Normalized Parameter = 1 − ( (Xi − Xref )2 σ2 X ) (3.1) where, Xi is the simulated result of the quantity (Head or Power) for the ith im- peller configuration with balancing holes at BEP (Q∗ = 1.0) condition, Xref is the simulated result of the same quantity for the reference impeller without balancing holes at BEP (Q∗ = 1.0) condition, and σX is the standard deviation of the sim- ulated results (Head or Power) across the 12 initial configurations generated using the Sobol sequence at BEP (Q∗ = 1.0) condition. To represent this issue more clearly, two extreme results from the BEP (Q∗ = 1.0) condition were considered. These results are presented in Table 3.5. The head (Hi) was normalized with the head (Href) for the BEP condition. Table 3.5: Illustration of normalization error using standard deviation. Configuration Hi Href (%) σX Normalized Data 1 100.03 0.129 0.9759 7 99.37 0.129 -7.2222 This example highlights how minor differences in head can yield disproportionately skewed normalized values. A similar behavior was observed for power. To address this, re-normalization of the already normalized data using min-max scaling was explored. This method maps the data into a [0, 1] range using the transformation: 41 3. Methodology x′ = x − min(x) max(x) − min(x) (3.2) Applying this to the normalized head values in Table 3.5, the re-normalized values become: x′ Config 1 = 0.9759 − (−7.2222) 0.9759 − (−7.2222) = 1.000, x′ Config 7 = −7.2222 − (−7.2222) 0.9759 − (−7.2222) = 0.000 Although this re-scaling constrains the data within the desired range, it does not cor- rect the core distortion. The original normalization already exaggerated differences due to the small standard deviation. As a result, min-max normalization merely rescales an already skewed distribution, compressing real variations and flattening intermediate configurations to values near 0 or 1. Therefore, this re-normalization method was ultimately not adopted, as it further reduced the interpretability of the normalized data and was inconsistent with the goal of constructing a smooth and meaningful target function. A more practical method was adopted: normalizing the data using a fixed percentage of the reference value (the impeller without balancing hole). The new normalization formula used is: Normalized Parameter = 1 − ( (Xi − Xref)2 (Xref · δ)2 ) (3.3) where, δ is the fixed percentage. δ = 0.2 % and δ = 1 % are used for power and head respectively. This approach ensured that the normalized data consistently fall within the range [0, 1], preserved the natural variation between samples, and enabled the construc- tion of a more robust and interpretable multi-objective target function. Hence, all subse