Slotting a hydrofoil is an effective way to passively control the cavitating flow to reduce the cavity pocket size leading to a reduction in vibration, noise, and erosion. However, this comes with losses in the hydrodynamic performance of the hydrofoil including its lift coefficient. To avoid this as much as possibl e, optimizing the slotted hydrofoil in terms of the location and the angle is of prime importance. The optimization is achieved by designing a deep Artificial Neural Network (ANN) to act as a surrogate model in the process of the genetic algorithm. The tra ining dataset in deep ANN is gathered through simulations that are based on a newly developed Frink finite volume solver for the preconditioned Euler equations. Results obtained indicate that this optimization approach is effective such that the cavity poc ket size can be shrieked by about 60 % with a penalty of about 10 % reduction in the hydrofoil’s lift coefficient. Keywords Optimized s lotted hydrofoil , Deep Artificial Neural Network (ANN) , Genetic algorithm, Surrogate model, Frink finite volume method , Cavitating flow 1 mahya. hajihassanpour@jacobs.com 2 1. Introduction Cavitating flows are widely seen in many practical applications such as in pumps, pipes, hydro turbines, and blood vessels, or on watercraft including ships, boats, amphibious, submarines, and hovercraft. Cavitation takes place when the liquid is accelerated and consequently the local pressure becomes lower than the vapor pressure. In such a condition, the preexisted vapor nucleus will grow in the liquid medium and eventually create a cavity pocket filled with both liquid and many vapor nucleus. When these nuclei arrive at places with higher pressures, condensation occurs. In the condensation process, the nu cleus violently collapses and high -pressure waves will be generated that cause erosion, annoying sounds, and vibration that may cause a shorter lifetime for the product [1]. Flow control techniques can be used to avoid these. Flow control techniques can be categorized into passive and active methods. In the case of active flow control, different methods such as the injection [2]–[5] or suction [6] are employed . As for the passive control methods vortex generators [7], [8] , leading -edge slats [9], leading -edge serration [10], slotted hydrofoils [11]–[15], micro -cylinders [16], hydrofoils with cavity [17], roughness [18], and flexible hydrofoil [19] have been used . Owing to the simplicity and economic efficiency of passive control methods over the active ones, former s are more applied to practical applications. Here, the passive flow control of the slotted hydrofoil is utilized . A few studies about the slotted hydrofoil have been performed in the literature [11]–[15]. The idea of using slotted hydrofoils comes from the work of Capurso et al. [11], [12] where three slots located near the leading -edge of the NACA 0009 hydrofoil have been considered . Their results show a 93% reduction for the total length of the cavity , but a 25% reduction in the lift coefficient. The optimization has not been performed by them. Conesa and Liem [13], [14] have studied the effect of the entry and exit widths of the slo t and its angle on the performance of the Clark -Y hydrofoil, but they have fixed the slot location near the lead ing-edge in the optimization process. Their study shows that the optimized slotted hydrofoil can achieve higher performance with a 50% reduction in cavity pocket in comparison with the base hydrofoil. An experimental and numerical study has been done by Ni et al. [15] to investigate the performance of a slotted hydrofoil operating close to a free surface. The slotted hydrofoi l is constructed from the base hydrofoil NACA 63 4-021. No optimization has been performed, and the slotted hydrofoil had a better performance in comparison with the base hydrofoil. Here , the optimization is done by changing 3 the location and the angle of the slot on the NACA 66 (MOD) base hydrofoil as it is widely used in many practical real-world artifacts . Nowadays, multi -objective optimizations are frequently used to design engineering artifact s. The solution to these complex design optimization problems, such as the problem that is addressed in this paper, are generally computationally expensive. To reduce the computational costs , surrogate model approximations are intensively used in t he literature [20], [21] . surrogate model s use limited combinations of inputs (sample points) and outputs to provide a mapping from the input space to the output space which essentially reduces the computational cost [22], [23] . surrogate model s come in several forms. Surrogate model s can be either deterministic or stochastic, low -order or high -order polynomial regression models, discrete or con tinuous, and local or global. Some of the better -known surrogate models are the Multivariate Adaptive Regression Splines, Radial Basis Function, Polynomial Regression, Neural Networks, Kriging, etc. [20]. The advantages or disadvantages of surrogate models vary c ase by case and we cannot easily say essentially one of them is better than the others, and an extensive study is needed between surrogate models to judge their priority, which is beyond the scope of this paper. In the current study, a deep Artificial Neur al Network (ANN) is employed . Choosing the sample training data points is of paramount importance to have a n accurate surrogate model that both represents the design space adequately and does not encounter any numerical difficulties. Usually, to determine the sample points the theory of the Design of Experiments (DoE) is employed. Several used DoE methods are the factorial, central decomposition, alphabetical optimal, and Box - Behnken methods [24], [25] . These methods usually tend to push the sample points towards the boundaries of the design space. More recent methods address this problem and try to uniformly distribute th e sample points within the boundaries. To name a few, Maximum Entropy, MiniMax, Orthogonal Arrays, and Latin Hypercube algorithms can be mentioned [20], [21], [24] . Latin Hypercube Sampling (LHS) is a straightforward wel l-established sampling method in the literature that is chosen to be used in this research. In LHS an even random sampling is employed and then all the chosen variables are randomly combined to give the samples required. While the surrogate models together with DoE methods are used to find out the optimization functions defined in the design space, the optimizers are required to figure out where minimums occur. The optimizers can be categorized into two gradient -based and gradient -free 4 groups [26]. As the name suggests, a gradient -based optimizer uses the gradient of the objective function to find the optimum solution. This dependence on the gradient raises two main issues. First, to use a gradient -based optimizer the gradient of t he function must be either known analytically or computed numerically. This increases the computational cost. Second, and more importantly, a gradient -based optimizer can get stuck in a local optimum. On the other hand, a gradient -free/heuristic/Metaheuris tics optimizer uses a random operator to search the design space. Decades of research resulted in the development of many Metaheuristics. To name a few well - known metaheuristics, we may mention Particle Swarm Optimization (PSO), Differential Evolution (DE) , Genetic Algorithms (GA), and Ant Colony Optimization (ACO) [27]. Inspired by Darwin’s natural evolution theory, the GA mimics the natural evolution to find the optimal solution [26]–[28]. Within the GA, the design parameters ar e interpreted as characteristics of individuals. Taking all individuals among a population, the solutions evolve through cross -over, mutation, and elitism operators to find the optimum. The random operator helps the Genetic Algorithm to explore the design space better , and the gradient -independence empowers it to deal with any kind of objective function. This rigorous searching capability grants the metaheuristics and the GA particularly the name “global optimization methods” [26] and it is used here. Realizations are needed to create a surrogate model (also called a surrogate model ). Cavitating flow realizations are usually achieved by simulating inviscid, laminar, or turbulent fluid flows. Assuming water as the main medium, the Reynolds number is often high in real applications so that the laminar cavitating flow can be only seen seldom, however, the inviscid flow condition can provide acceptable results and the turbulent flow simulations are usually expensive . The pressure distribution over the hydrofoil determines the cavity pocket size and the lift coefficient can be reasonably captured by the solution of inviscid Euler equations. In an inviscid flow simulation, a fewer number of grid points in comparison with turbulent flows are required in simulations and then the computational cost can be dramatically decreased. Therefore, the preconditioned Euler equations representing inviscid flows are used here. The preconditioned Euler equations should be solved using an appropriate numerical method. Different methods have been utilized in the literature for solving these governing equations, such as the central finite difference [29], compact [30], central finite volume [31], Quadratic Upstream Interpolation for the Convective Kinematics (QUICK) [32], weighted 5 essentially non -oscillating (WENO) [33], immersed boundary [34], upstream finite element [35], and discontinuous Galerkin [36]–[38] methods. Here, the Frink numerical method [39] is applied to discretize the governing equation, and an appropriate artificial dissipation term, based on the one proposed by Jameson et al. [40], is added f or suppressing artificial oscillations in the solution domain and to provide a stable , accurate , and reasonably fast solver. In this paper, the algorithm of the optimizer is given in Section 2. The design of the deep ANN is described in Section 3 and the details of the newly Frink finite volume -based simulator is given in Section 4. Results obtain are discussed in Section 5 and the paper is ended up through a conclusion in Section 6. 2. Optimizer The optimization parameters are the slot angle 𝜃 and its location 𝑥𝑠𝑙𝑜𝑡 with constraints of 0°≤ 𝜃≤85° and 0≤𝑥𝑠𝑙𝑜𝑡≤0.4. The optimization parameters are geometrically introduced in Fig. 1. As shown , the slot location varies from the middle chord up to the leading -edge and its angle is defined clockwise. The Latin Hypercube Sampling method (LHS) with an affine mapping is used to generate initial sampling points 𝑁𝑠𝑎𝑚𝑝𝑙𝑒𝑠 in this (𝜃, 𝑥𝑠𝑙𝑜𝑡) design space. The initial number of sampling points are supposed to be 4 0 points. Then, the CFD realizations (Section 4) are performed for the given sampling points to calculate the lift coefficient and the total length of the cavity appeared on the hydrofoil . Therefore, the lift coefficient and the total leng th of the cavity are our two objective function s. The lift coefficient function is represented by 𝐶𝑙=𝑓(𝑥𝑠𝑙𝑜𝑡,𝜃) and the total length of cavity is denoted by 𝐿𝑐=𝑔(𝑥𝑠𝑙𝑜𝑡,𝜃) will be found using a deep ANN (Section 3) which is exploited to construct the surrogate model . The lift coefficient can be calculated using the integral of the pressure in the direction perpendicular to the far -field flow divided by 0.5𝜌𝑢∞2. In the case of the total length of the cavity 𝐿𝑐, it will be calculated using the following relation : 𝐿𝑐=𝐿𝑐1+𝐿𝑐2 (1) where 𝐿𝑐1 denotes the length of the leading -edge cavity and 𝐿𝑐2 is the length of the mid -chord cavity. The leading -edge cavity is a cavity formed on the front hydrofoil, and the mid -chord cavity is a cavity that appeared on the rear hydrofoil created by slotting the original hydrofoil. The length of the leading -edge/mid -chord cavity is the horizontal distance between two points where the 6 pressure over the hydrofoil becomes lower than the vapor pressure and where it becomes higher than the vapor pressure. The two objective function s of 𝐶𝑙=𝑓(𝑥𝑠𝑙𝑜𝑡,𝜃) and 𝐿𝑐=𝑔(𝑥𝑠𝑙𝑜𝑡,𝜃) will serve as an input to a genetic algorithm (GA) to result in an accurate and fast optimization scheme. The Pareto front comes out from GA with 𝑁𝑝𝑎𝑟𝑡𝑜 design point s. These procedures are repeated 𝑁𝐺𝐴 times. In each iteration, the 𝑁𝑝𝑎𝑟𝑡𝑜 new design points will be added to the sample points , i.e., 𝑁𝑠𝑎𝑚𝑝𝑙𝑒𝑠 will be updated to 𝑁𝑠𝑎𝑚𝑝𝑙𝑒𝑠+𝑁𝑝𝑎𝑟𝑡𝑜. Subsequently, CFD realizations , deep ANN, and GA will be carried out. Here, 𝑁𝐺𝐴=5 is set which means that five generations are considered in the genetic algorithm. The number of new design points will not be greater than 𝑁𝑝𝑎𝑟𝑡𝑜,𝑚𝑎𝑥=20 in each iteration, and if not, 20 points will be selected randomly. The latter is presumed t o reduce the computational cost. This optimization procedure/algorithm is illustrated in Fig. 2. This optimization procedure is performed in Matlab. Th e information of the sampling points will be automatically transformed into a meshing tool which is Gambit to generate the mesh for the CFD tool. The CFD tool is also a homemade code using the Frick method and C++ as a programming language (Section 4). Fig. 1 Nomenclatures related to the optimization parameters, the total length of the cavity, the length of the leading -edge cavity, and the length of the mid -chord cavity. xy -0.4 -0.2 0 0.2 0.4-0.4-0.200.20.4 xslot Lc2 Lc1Lc=Lc1+Lc2 7 Fig. 2 Optimizer f lowchart. 8 3. Surrogate model (deep ANN) Knowing 𝑁𝐿𝐻𝑆 input values for 𝜃 and 𝑥𝑠𝑙𝑜𝑡 using LHS method and their association target output values of 𝐶𝑙=𝑓(𝑥𝑠𝑙𝑜𝑡,𝜃) and 𝐿𝑐=𝑔(𝑥𝑠𝑙𝑜𝑡,𝜃) using the CFD tool , a deep ANN can be trained . The training dataset is created by 80% of data while the remaining 20% of data is used to verify the trained network. The preprocessing on the training dataset should be done before training a deep ANN. This preprocessing includes the normalization of data to the interval of [0,1]. Normalization increases the accuracy and speed of the ANN. The first layer of the designed ANN is an input layer with two features (𝜃 and 𝑥𝑠𝑙𝑜𝑡,). This layer is connected to a triple -layer that will be repeated 𝑁𝑟 times . The first layer in this triple layer is a fully connected layer with 𝑁𝑛 neurons. The next one is the leaky Rectified Linear Unit (ReLU) that can be efficiently employed for the regression jobs. The last one in the triple -layer is the batch normalization layer aiming to reduce errors of the ANN and increase its accuracy by retaining numbers throughout the ANN on the same scale. After repeating this triple -layer 𝑁𝑟 times, then a fully connected layer is used to give out the data to the last layer which is the regression layer to estimate 𝐶𝑙=𝑓(𝑥𝑠𝑙𝑜𝑡,𝜃) and 𝐿𝑐=𝑔(𝑥𝑠𝑙𝑜𝑡,𝜃). A schematic of the designed deep ANN is shown in Fig. 3. Fig. 3 Architecture of the designed ANN network A residual minimizer is performed to properly adjust the weights and biases of the deep ANN . The residual minimizer tries to minimize a loss function indicating the differences between the deep ANN outputs (estimated surrogate models at the sample points) and the target outputs (CFD results) . The Adam method is selected as the residual minimizer . The number of epochs in the Adam method is 250 and the learning rate equals 0.08. Our study shows that thes e values for epochs and learning rate leads to an accurat e and fast deep ANN. The number of repeated triple -layers is 𝑁𝑟=4 and the number of neurons is 𝑁𝑛=30, which plus prescribed deep ANN settings leads to not getting stuck on overfitting and underfitting problems. 9 4. Simulator Steady inviscid cavitating fluid flows can be mathematically described using preconditioned Euler equations with the artificial compressibility method as follows: 𝚷𝜕𝑸 𝜕𝜏+𝜕𝑭 𝜕𝑥+𝜕𝑮 𝜕𝑦−𝑺𝒄=0 (2) where 𝚷 indicates the preconditioning matrix which is 𝚷= [ 1 𝜌𝑚𝛽2000 0𝜌𝑚0𝑢Δ𝜌𝑙 00𝜌𝑚𝑣Δ𝜌𝑙 𝛼𝑙 𝜌𝑚𝛽2001] (3) The artificial time is denoted by 𝜏 and derivatives with respect to this variable should go to zero for satisfying the consistency and convergence conditions. This means that this formulation can be used only for steady -state solutions. This is corresponding with the physics we are dealing with because we are going to optimize a hydrofoil experiencing quasi -steady sheet cavitation flow. The cavitation source term 𝑺𝒄 plays an important rule in cavitation modeling and it is defined in a way to appropriately model the mass transfer between the liquid and vapor phases that occurred in the cavitation process. We will further discuss this term in the next section. The flux vectors 𝑯= (𝑭,𝑮) and the solution vector 𝑸 are written as follows: 𝑸=[𝑝 𝑢 𝑣 𝛼𝑙], 𝑭=[𝑢 𝜌𝑚𝑢2+𝑝 𝜌𝑚𝑢𝑣 𝛼𝑙𝑢], 𝑮=[𝑣 𝜌𝑚𝑢𝑣 𝜌𝑚𝑣2+𝑝 𝛼𝑙𝑣] (4) The solution vector is consisting of the pressure 𝑝, the velocity vector (𝑢,𝑣), and the liquid volume fraction 𝛼𝑙. The subscribes 𝑙, 𝑣, and 𝑚 are added to facilitate identifying the liquid, vapor, and mixture states. While the liquid density 𝜌𝑙 and 𝜌𝑣 are supposed to be constant, the mixture density is allowed to vary according to the relation 𝜌𝑚=𝛼𝑙𝜌𝑙+𝛼𝑣𝜌𝑣. Note that the summation of volume fraction of different phases is always unity, i.e., 𝛼𝑙+𝛼𝑣=1. Moreover, 𝜌𝑙=1 and 𝜌𝑣=0.01 are used here that prove to be suitable values for simulating cavitating flows in a water medium. In Eq. (3), the density jump is indicated by Δ𝜌𝑙=𝜌𝑙−𝜌𝑣. 10 Following the concept of the artificial compressibility method, the coefficient denoted by 𝛽 is used in the governing equations to tackle the pressure -velocity coupling problem pertaining to incompressible flows. It is a user -defined parameter and here we simply use a constant value 𝛽= 3 similar to [36], [41] . 4.1. Cavitation modeling As mentioned in the previous section, the cavitation process is modeled by the cavitation source term 𝑆𝑐 which is 𝑺𝒄= [ (𝑚̇++𝑚̇−)(1−1 𝜌𝑣) 0 0 𝑚̇++𝑚̇−] (5) where 𝑚̇− models the evaporation phenomena which is the mass transfer from the liquid to the vapor phase, and it is vice versa for 𝑚̇+ which represents the condensation process wherein the vapor phase transits to the liquid phase. These two variables are differently defined in the literature [42]–[46]. Here, the relations proposed by Merkel et al. [42] is exploited, i.e., 𝑚̇+=(𝐶𝑝𝑟𝑜𝑑 𝑡∞)(1−𝛼𝑙)𝑀𝑎𝑥(0,𝑝−𝑝𝑣) 𝑚̇−=(𝐶𝑑𝑒𝑠𝑡 𝑡∞)1 𝜌𝑣 𝛼𝑙𝑀𝑖𝑛(0,𝑝−𝑝𝑣) (6) These relations have been proposed by Merkel et al. [42] that derived from their experimental work. N umerical studies [36], [47] show that these relations are reliable and accurate for sim ulating cavitating flows. The parameters 𝐶𝑝𝑟𝑜𝑑/𝑡∞ and 𝐶𝑑𝑒𝑠𝑡/𝑡∞ are the constant parameters which control the liquid production and destruction rates, respectively. They are user -defined parameters where 𝐶𝑝𝑟𝑜𝑑/𝑡∞=1 and 𝐶𝑑𝑒𝑠𝑡/𝑡∞=80 are used here. Numerical studies done in literature [36], [47] prove these values lead to accurate numerical results. Finally, the vapor pressure 𝑝𝑣 is obtained according to the predefined cavitation number 𝜎=2(𝑝∞−𝑝𝑣) where 𝑝∞ is the pressure far enough from the hydrofoil surface whic h is assumed unit. Note that all the variables used in the present study are non -dimensional identical to the ones presented in [36]. 11 4.2. Spatial discretization So that being able to simulate cavitating flows, we should suitably address the solution of the governing equation (2) by means of a numerical method. H ere, the Frink numerical method is used to discretize the spatial derivative terms in the governing equations. By integrating Eq. (2) through the surface A of an element, it can be rewritten as ∫[𝚷𝜕𝑸 𝜕𝜏+𝜕(𝑭) 𝜕𝑥+𝜕(𝑮) 𝜕𝑦−𝑺𝒄]𝑑A 𝐀=0 (7) In the finite volume method variables are assumed to be constant through an element, then ∫𝚷𝜕𝑸 𝜕𝜏𝑑A 𝐀=𝚷𝜕𝑸 𝜕𝜏A (8) For the spatial derivative of inviscid fluxes, we have ∫[𝜕𝑭 𝜕𝑥+𝜕𝑮 𝜕𝑦]𝑑A 𝐀=∫[𝛁∙𝑯]𝑑A 𝐀 (9) and then using the divergence theorem applied to triangular domain results in ∫[𝛁.𝑯]𝑑A 𝐀=∮[𝑯∗.𝒏̂] 𝑑𝐿 𝛛𝐀=∑(𝑯∗.𝒏̂)𝒌𝐿𝑘𝑵𝒇𝒂𝒄𝒆𝒔=𝟑 𝒌=𝟏 (10) where 𝑯𝒌∗ is the numerical flux for the 𝑘-th face of the element, and it depends on the flux vectors in two neighboring elements sharing that face. After finding the flux vectors in two neighboring elements by exploiting the Frink numerical method, the Lax numerical flux method is used to calculate the numerical flux. The calculation of the Lax numerical flux me thod will be explained later. Variable 𝒏̂ is the normal outward vector. The face length of the triangle element is denoted by 𝐿. The integration of the cavitation source term 𝑺𝒄 can be written as ∫𝑺𝒄 𝑑A 𝐀=𝑺𝒄 A (11) 12 The integration of the artificial dissipation terms should be carefully evaluated, and it will be addressed in the next section. Finally, the discretization process, by substituting Eqs. (8)-(11) into (7), results in 𝜕𝑸 𝜕𝜏=𝚷−𝟏 𝐴[𝑺𝒄 A−∑(𝑯∗.𝒏̂)𝒌𝐿𝑘𝑵𝒇𝒂𝒄𝒆𝒔 𝒌=𝟏]=𝑹 (12) where 𝑅 is called the right -hand side vector. 4.2.1. The Frink numerical method The higher -order finite volume methods can be derived using reconstructing fluxes at the element faces. To show how the Frink numerical method can be used to increase the order of accuracy of the solution from the first -order in the classical finite volume method to the second -order one, refer to Fig. 4. While in the classical first -order finite volume methods variables at the faces can be read as 𝑞𝑓1=𝑞𝑓2=𝑞𝑓3=𝑞𝑐, in the Frink numerical method [39], they are reconstructed as follows: 𝑞𝑓1=𝑞𝑐+1 3[1 2(𝑞𝑛1+𝑞𝑛2)−𝑞𝑛3] 𝑞𝑓2=𝑞𝑐+1 3[1 2(𝑞𝑛2+𝑞𝑛3)−𝑞𝑛1] 𝑞𝑓3=𝑞𝑐+1 3[1 2(𝑞𝑛1+𝑞𝑛3)−𝑞𝑛2] (13) The remaining variables to be defined are the nodal values denoted by 𝑞𝑛1, 𝑞𝑛2, and 𝑞𝑛3. An area - weighted average, among the elements sharing the node, is used to calculate values in nodes. 13 Fig. 4 Nomenclatures used for constructing fluxes in the Frink numerical method. 4.2.2. Numerical flux The Lax numerical flux is used here to address the Riemann problem for the jumps between two neighboring elements. If the flux in a face is denoted by the subscript 𝐿 and its neighbor’s value by 𝑅, then 𝑯𝒌∗=𝟏 𝟐(𝑯𝑳+𝑯𝑹−𝐶𝜆𝚷(𝑸𝑹−𝑸𝑳)) (14) where 𝜆 is the largest eigenvalue of the system of equations 𝜆=√𝑢𝐿2+𝑣𝐿2+√𝑢𝐿2+𝑣𝐿2+𝛽2 (15) and 𝐶 is a constant parameter with 𝐶≤1 for satisfying stability conditions. Selecting higher values for 𝐶 means adding more dissipation to the solution domain that may cause reducing global accuracy. Here, 𝐶=0.25 is used to consequently provides accuracy and stability needs. qf3qcqn3qf2 qf1 qn1qn2 14 4.2.3. Artificial dissipation Although discretization has been done, we should care about the artificial non -physical oscillations that occur near the cavity zone which could make the solution uns table. To prevent such undesirable oscillations in the solution, the artificial dissipation term 𝑫 [40] is added to the right - hand side of Eq. (12) as follows: 𝜕𝑸 𝜕𝜏=𝑹+𝑫=𝑹′ (16) The artificial dissipation for the 𝑗-th element can be written as: 𝑫𝒋=∑𝒅𝑗𝑘𝑁𝑓𝑎𝑐𝑒𝑠=3 𝑘=1 (17) where 𝑘 indicates the face number. If the solution vector in neighbor element through the 𝑘-th face is shown by 𝑄𝑗𝑘, then the dissipation term can be written as follows: 𝒅𝒋𝒌=−(𝐴𝑗 ∆𝜏+𝐴𝑗𝑘 ∆𝜏)[𝜀𝑗𝑘 2(𝑸𝒋−𝑸𝒋𝒌)] (18) where ∆𝜏 is the artificial time step. The artificial viscosity coefficient 𝜀 has to be determined such that tends to zero where the solution is smooth and to a definite value in non -smooth regions which normally occur near the cavity zone which is the origin of the artifact oscillations generated due to the sharp jump in the density variable across the interface of the two phases. In order to determine the artificial viscosity coefficient, we follow the idea presented in Ja meson et al. work [40] which is: 𝜀𝑗𝑘=𝜅𝜌𝛾𝑗𝑘 𝛾𝑗𝑘=|𝜌𝑗−𝜌𝑗𝑘 𝜌𝑗+𝜌𝑗𝑘| (19) and 𝜅𝜌=0.025 is chosen here. 4.3. Temporal discretization The Euler method is used here for discretizing the temporal derivative in Eq. (16) 15 𝑄𝑚+1−𝑄𝑚 ∆𝜏=𝑹′𝒎 (20) where superscript 𝑚 indicates the time -level. 5. Results In this section, the numerical results are presented. The simulator is first evaluated in subsection 5.1 wherein the sheet /leading -edge and mid -chord cavitations have been simulated and are compared with the experimental data. After proving that the numerical method works fine, then the results for the slotted hydrofoil and its optimized shape will be discussed. 5.1. Validation Once the original hydrofoil NACA 66(MOD) is slotted, instead of one cavity pocket over the original hydrofoil surface, two -cavity pockets might form over the modified hydrofoil. The first cavity pocket over the modified hydrofoil can be seen as a sheet /leading-edge cavity , and the second one is like the mid -chord cavity over the original hydrofoil. Thus, the sheet and mid -chord cavitation phenomenon are simulated here to show that the numerical method arranged here can be effectively used for the optimization of the slotted hydrofoil. The sheet cavitation will be created over the original hydrofoil when water flows over NACA 66 (MOD) at the angle of attack 𝛼=4° and cavitation number 𝜎=0.84, and by considering 𝛼=1° and 𝜎=0.43 the mid -chord cavitation will be formed. In Fig. 5, the generated grid with 16330 triangle elements is shown. The grid type is chosen to be triangular because they can be used over arbitrary and complex geometries such as the cases we are deal ing with in this paper, i.e., slotted hydrofoil. Moreover, they are suitably working in inviscid flows, because there is no need for having high -resolution grids close to the walls due to the absence of the boundary layer. As can be seen in Fig. 5, the radius of the far -field boundary is assumed ten times of the hydrofoil chord. Moreover, a finer grid is used close to the hydrofoil surface while a coarser grid is used far from its surface. Using high -resolution grids near the hydrofoil will help to sharply capture the cavity pocket, and a coarse grid near the far -field boundary will aid to damp the artificial incoming and outgoing waves occurred there. 16 In Fig. 6, the numerical results for the simulation of the sheet cavitation over the original NACA 66 (MOD) hydrofoil are illustrated. The density, pressure, and velocity magnitude contours indicate tha t the cavitation phenomenon can be expected over the suction side of the hydrofoil wherein the velocity magnitude is increased and subsequently the pressure is decreased to a lower value than the vapor pressure. Moreover, the velocity magnitude immediately reduces when the cavity pocket is closed which is a sign for the re -entrant jet. In Fig. 7, the calculated surface pressure coefficient by the present numerical method is compared with the experimental data [48]. As can be seen in this figure, the agreement between numerical results is satisfactory. The pressure in the cavity pocket remains unchangeable that excellently follows experimental pieces of evidence. In addition, the pressure recovery at the trailing edge is captured by the present numerical method efficiently. In Fig. 8 the density, pressure and velocity magnitude contours for the mid -chord cavitation type has been plotted. As shown the cavitation does not start from the lead ing-edge , but starts from the middle of the suction side. I n Fig. 9, the numerical results by utilizing the Frink numerical method are shown and compared with the experimental data [48] for this cavitation type . As shown in this figure, the numerical results obey the expected physical behaviors and also agree with the experimental data. From the obtained numerical results for the simulation of these two cavitation types, it can be concluded that the numerical framework set up here can be used to effectively simulate different quasi -steady cavitation flows encountered in the present study. 17 Fig. 5 Generated grid for the original NACA 66 (MOD) hydrofoil . xy -10 -5 0 5 10-10-50510 18 (a) (b) (c) Fig. 6 Contours of a) the density , b) the pressure, and c) the velocity magnitude with the streamlines computed by the Frink numerical method for the problem of cavitating flow over NACA 66 (MOD) at 𝜎=0.84 and 𝛼=4°. xy -0.5 0 0.5-0.200.20.40.6  0.94 0.84 0.74 0.64 0.54 0.44 0.34 0.24 0.14 xy -0.5 0 0.5-0.200.20.40.6 p 1.37 1.27 1.17 1.07 0.97 0.87 0.77 0.67 0.57 xy -0.5 0 0.5-0.200.20.40.6 Vel.Mag. 1.34 1.24 1.14 1.04 0.94 0.84 0.74 0.64 0.54 0.44 0.34 0.24 0.14 19 Fig. 7 Comparison of the pressure coefficient obtained by the Frink numerical method and experimental data [48] for the problem of cavitating flow over NACA 66 (MOD) at 𝜎=0.84 and 𝛼=4°. xCp -0.4 -0.2 0 0.2 0.4-1 -0.5 0 0.5 1Frink Experiment 20 (a) (b) (c) Fig. 8 Contours of a) the density, b) the pressure, and c) the velocity magnitude with the streamlines computed by the Frink numerical method for the problem of cavitating flow over NACA 66 (MOD) at 𝜎=0.43 and 𝛼=1°. xy -0.5 0 0.5-0.200.20.40.6  0.99 0.89 0.79 0.69 0.59 0.49 0.39 xy -0.5 0 0.5-0.200.20.40.6 p 1.38 1.28 1.18 1.08 0.98 0.88 0.78 xy -0.5 0 0.5-0.200.20.40.6 Vel.Mag. 1.12 1.02 0.92 0.82 0.72 0.62 0.52 0.42 0.32 0.22 0.12 21 Fig. 9 Comparison of the pressure coefficient obtained by the Frink numerical method and experimental data [48] for the problem of cavitating flow over NACA 66 (MOD) at 𝜎=0.43 and 𝛼=1°. 5.2. Optimization of slotted hydrofoil The goal of this study is to optimize the slotted NACA 66 (MOD) hydrofoil immersed in water at cavitation number 𝜎=0.84 and angle of attack 𝛼=4°. As discussed in detail in the optimization procedure, 40 initial sampling points are used to start the optimization process . These initial sampling points are obtained by employing the LHS method , and these points are depicted in Fig. 10(a). By advancing in the optimization process , new design points will be added to the initial sampling points to create and assess the objective functions and to finally find out the optimum points. All the points used in the optimization process are illustrated in Fig. 10(b) and it can be seen that many new points are inserted in the upper right side of the figure where the slots a re created close to the leading -edge with higher angles. It means that the optimizer seeks the optimal point there, and by adding many new points in that region, it improves the hig h-fidelity model. In Fig. 11, the contours related to 𝐿𝑐1 and 𝐿𝑐2 obtained using the deep ANN are shown. Form Fig. 11(a) it can be concluded that generally at a given slot location the leading -edge cavity xCp -0.4 -0.2 0 0.2 0.4-0.4 -0.2 0 0.2 0.4 0.6Frink Experiment 22 will reach its minimum size around 𝜃=45°, while lower an d higher angles than this angle w ill cause having larger leading -edge cavity size. At a given slot angle, moving the slot location from the middle to the leading -edge will favorably result in decreasing in 𝐿𝑐1. Thus, smaller 𝐿𝑐1 will be found around 𝜃=45° and high slot angles. In the case of 𝐿𝑐2, as depicted in Fig. 11(b), at a given slot location the mid -chord cavity will r each its maximum size around 𝜃=60°, while lower and higher angles than this angle w ill cause having larger leading -edge cavity size. At a given slot angle, by increasing the slot location , 𝐿𝑐2 will disadvantageously increase. Thus, for having smaller 𝐿𝑐2, higher slot angles should be used with a slot placed around the middle of the hydrofoil . The contours of the objective functions 𝐶𝑙 and 𝐿𝑐 are illustrated in Fig. 12. The objective functions are purely non -linear with respect to optimization parameters. This figure can be very useful for engineers when they are designing artifacts to find out a design point passing their presumed criteria and to figure out how changing the design point would affect the performance of the artifact. According to the results plotted in Fig. 12, we can roughly argue that there are two domains where 𝐿𝑐 becomes minimum 𝐷𝑜𝑚(𝐿𝑐𝑚𝑖𝑛)1={(𝜃,𝑥𝑠𝑙𝑜𝑡)∈ℝ|25°<𝜃<65° ,0<𝑥𝑠𝑙𝑜𝑡<0.15} 𝐷𝑜𝑚(𝐿𝑐𝑚𝑖𝑛)2={(𝜃,𝑥𝑠𝑙𝑜𝑡)∈ℝ|72.5°<𝜃<85° ,0.1<𝑥𝑠𝑙𝑜𝑡<0.4} (21) and there is one domain that 𝐶𝑙 becomes maximum 𝐷𝑜𝑚(𝐶𝑙𝑚𝑎𝑥)={(𝜃,𝑥𝑠𝑙𝑜𝑡)∈ℝ|40°<𝜃<85° ,0.15<𝑥𝑠𝑙𝑜𝑡<0.4} (22) Then, it is expected that the optimized points placed in the set stated hereunder: (𝐷𝑜𝑚(𝐶𝑙𝑚𝑎𝑥)∩𝐷𝑜𝑚(𝐿𝑐𝑚𝑖𝑛)1)∪(𝐷𝑜𝑚(𝐶𝑙𝑚𝑎𝑥)∩𝐷𝑜𝑚(𝐿𝑐𝑚𝑖𝑛)2) ={(𝜃,𝑥𝑠𝑙𝑜𝑡)∈ℝ|72.5°<𝜃<85° ,0.15<𝑥𝑠𝑙𝑜𝑡<0.4} (23) which is exactly corresponding with the domain that the optimizer seeks for the optimal points as depicted in the right upper part of Fig. 10(b). It is evident that the optimization process is suitably designed and works appropriately. Before presenting the optimal points, it is important to assess the accuracy of the surrogate mode l and its validity. For this aim, at different points than the trained dataset, along the midlines of the design space, the CFD results are compared with the data 23 extracted from the surrogate model that are given in Table 1. The percentage of the error is defined as 𝐸𝑟𝑟𝑜𝑟(𝐿𝑐)=|𝐿𝑐𝐶𝐹𝐷−𝐿𝑐𝑠𝑢𝑟𝑟𝑜𝑔𝑎𝑡𝑒| 𝐿𝑐𝐶𝐹𝐷×100 𝐸𝑟𝑟𝑜𝑟(𝐶𝑙)=|𝐶𝑙𝐶𝐹𝐷−𝐶𝑙𝑠𝑢𝑟𝑟𝑜𝑔𝑎𝑡𝑒| |𝐶𝑙𝐶𝐹𝐷|×100 (24) and it is seen that the maximum percentage of the error in predicting 𝐿𝑐 is lower than 4.06 percent, while for 𝐶𝑙 this error is lower than 2.24. These amount s of errors are in an acceptable range, and it can be concluded that the surrogate model is trained properly and is reliable. In Figs. 13 and 14, the density and pressure contours for points used in Table 1 are plotted. From Fig. 13 and Table 1 it can be deduced that for 𝜃=42.5° and by pushing the slo t from the middle chord to the leading - edge, the lift coefficient is desirably increased while the total length of cavity is unfavorably increased. However, form Fig. 14 and Table 1 it can be seen that for 𝑥𝑠𝑙𝑜𝑡=0.2 and by increasing the slot angle, while 𝐿𝑐 is always desirably decreasing, 𝐶𝑙 is first decreasing and then increasing. In Table 2, the output of t he optimization process , which is the final Pareto front , is given. While any point in the Pareto front is a candidate to be served as an optimal point, another criterion should be added to help the optimal point selection. This criterion is arbitrary. Her e, we are seeking a point that shows the best performance in comparison with the average values of all points in the Pareto front. These average values are indicated by 𝐿𝑐̅̅̅ and 𝐶𝑙̅ and are presented in the last row. In other words, the optimal point would be a point with max (𝐿𝑐̅̅̅−𝐿𝑐) and max (𝐶𝑙−𝐶𝑙̅) which is achieved for (𝜃,𝑥𝑠𝑙𝑜𝑡)=(80.461,0.354). In Fig. 15, generated gird for this optimal slotted hydrofoil is shown where, similar to the original hydrofoil illustrated in Fig. 5, a finder grid resolution is utilized close to the hydrofoil surface in comparison with the one s far enough from it. The pressure contour with streamlines and the density contour is shown in Fig. 16. It can be seen that by the slot the fluid with higher -pressure on the lower surface of the hydrofoil is transferred to the suction side with lower pressure resulting in shrinking the cavity pocket. In fact, the high -pressure fluid is injected to the suction side with lower -pressure and prevents the pressure becoming lower than the vapor pressure and subsequently results in not to happen cavitation. However, injecting a high -pressure fluid over the 24 hydrofoil will decrease the lift coefficient. This phenomenon can be clearly seen in Fig. 17 where the pressure coefficient over the surface of the optimized slotted hydrofoil is compared with the original hydrofoil. In this way, we can favorably shrink the cavity pocket, but the lift coefficient will be disadvantageously decreased. Basically, imposing higher pressure over the upper surface of the hydrofoil will cause a downward force that means decreasing in the lift. The total length o f the cavity that appeared on the suction side for the optimized slotted hydrofoil is 𝐿𝑐=0.2006 , while it is 0.482 for the original hydrofoil that shows an improvement. In the case of the lift coefficient, it equals 0.6098 for the optimized slotted hydro foil while it is 0.688 for the original hydrofoil that shows a deterioration. These result s indicate that with a penalty of about 11.4% in the lift coefficient, it is possible to decrease the cavity size by about 58.4%. 25 (a) (b) Fig. 10 a) Initial sampling using 40 points and b) final training samples. xslot 0 20 40 60 8000.050.10.150.20.250.30.350.4 xslot 0 20 40 60 8000.050.10.150.20.250.30.350.4 26 (a) (b) Fig. 11 Contours for a) 𝐿𝑐1 and b) 𝐿𝑐2 obtained using the fitting method . xslot 0 20 40 60 8000.050.10.150.20.250.30.350.4 Lc1 0.19 0.17 0.15 0.13 0.11 0.09 0.07 0.05 0.03 0.01 xslot 0 20 40 60 8000.050.10.150.20.250.30.350.4 Lc2 0.36 0.32 0.28 0.24 0.20 0.16 0.12 0.08 0.04 0.01 27 (a) (b) Fig. 12 Contours for a) 𝐿𝑐 and b) 𝐶𝑙 obtained using the fitting method. xslot 0 20 40 60 8000.050.10.150.20.250.30.350.4 Lc 0.40 0.37 0.34 0.32 0.29 0.27 0.24 0.22 0.20 0.17 0.14 0.12 0.10 xslot 0 20 40 60 8000.050.10.150.20.250.30.350.4 Cl 0.62 0.61 0.60 0.59 0.58 0.57 0.56 0.55 0.54 0.53 0.52 0.51 0.50 0.49 0.48 28 Table 1 Investigation of the accuracy of the surrogate model 𝜃 𝑥𝑠𝑙𝑜𝑡 CFD Surrogate model Error (%) 𝐿𝑐 𝐶𝑙 𝐿𝑐 𝐶𝑙 𝐿𝑐 𝐶𝑙 42.5° 0.00 0.118 0.476 0.123 0.484 4.778 1.591 42.5° 0.20 0.231 0.519 0.222 0.514 3.723 0.966 42.5° 0.40 0.323 0.564 0.339 0.577 4.978 2.244 0.0° 0.20 0.284 0.574 0.297 0.575 4.532 0.068 42.5° 0.20 0.231 0.519 0.222 0.514 3.723 0.966 85.0° 0.20 0.157 0.584 0.162 0.586 2.852 0.286 29 (a) (b) (c) Fig. 13 The pressure (left column) and density (right column) contours for 𝜃=42.5° and a) 𝑥𝑠𝑙𝑜𝑡=0.0, b) 0.2, and c) 0.4 xy -0.5 0 0.5-0.200.20.40.6 p 1.38 1.28 1.18 1.08 0.98 0.88 0.78 0.68 0.58 xy -0.5 0 0.5-0.200.20.40.6  0.94 0.84 0.74 0.64 0.54 0.44 0.34 0.24 0.14 xy -0.5 0 0.5-0.200.20.40.6 p 1.38 1.28 1.18 1.08 0.98 0.88 0.78 0.68 0.58 xy -0.5 0 0.5-0.200.20.40.6  0.94 0.84 0.74 0.64 0.54 0.44 0.34 0.24 0.14 xy -0.5 0 0.5-0.200.20.40.6 p 1.38 1.28 1.18 1.08 0.98 0.88 0.78 0.68 0.58 xy -0.5 0 0.5-0.200.20.40.6  0.94 0.84 0.74 0.64 0.54 0.44 0.34 0.24 0.14 30 (a) (b) (c) Fig. 14 The pressure (left column) and density (right column) contours for 𝑥𝑠𝑙𝑜𝑡=0.2 and a ) 𝜃=0.0°, b) 42.5°, and c ) 85.0°. xy -0.5 0 0.5-0.200.20.40.6 p 1.38 1.28 1.18 1.08 0.98 0.88 0.78 0.68 0.58 xy -0.5 0 0.5-0.200.20.40.6  0.94 0.84 0.74 0.64 0.54 0.44 0.34 0.24 0.14 31 Table 2 Pareto front design points (GA output) 𝜃 𝑥𝑠𝑙𝑜𝑡 𝐿𝑐 𝐶𝑙 𝐿𝑐̅̅̅−𝐿𝑐 𝐶𝑙−𝐶𝑙̅ 84.545 0.003 0.1293 0.5976 0.0814 -0.0119 80.780 0.094 0.1945 0.6091 0.0162 -0.0005 79.936 0.192 0.2107 0.6112 0.0001 0.0017 78.943 0.138 0.2295 0.6132 -0.0187 0.0036 75.576 0.146 0.2882 0.6169 -0.0775 0.0073 80.461 0.354 0.2006 0.6098 0.0102 0.0003 Selected 79.473 0.284 0.2196 0.6123 -0.0089 0.0028 76.866 0.339 0.2676 0.6164 -0.0569 0.0068 82.397 0.111 0.1648 0.6046 0.0460 -0.0050 78.373 0.198 0.2402 0.6143 -0.0295 0.0047 77.263 0.247 0.2605 0.6158 -0.0497 0.0063 82.786 0.003 0.1578 0.6033 0.0529 -0.0063 82.999 0.047 0.1542 0.6027 0.0566 -0.0069 83.464 0.202 0.1465 0.6014 0.0643 -0.0082 78.811 0.302 0.2322 0.6136 -0.0215 0.0041 76.170 0.073 0.2786 0.6165 -0.0679 0.0070 84.996 0.001 0.1228 0.5961 0.0880 -0.0135 75.117 0.400 0.2958 0.6173 -0.0850 0.0077 𝐿𝑐̅̅̅= 0.2107 𝐶𝑙̅= 0.6096 32 Fig. 15 Generated grid for the optimal slotted NACA 66 (MOD) hydrofoil. xy -10 -5 0 5 10-10-50510 33 (a) (b) Fig. 16 The pressure (left column) and density (right column) contours for a) the original hydrofoil and for b) the optimized slotted hydrofoil . xy -0.5 0 0.5-0.200.20.40.6 p 1.37 1.27 1.17 1.07 0.97 0.87 0.77 0.67 0.57 xy -0.5 0 0.5-0.200.20.40.6  0.94 0.84 0.74 0.64 0.54 0.44 0.34 0.24 0.14 xy -0.5 0 0.5-0.200.20.40.6 p 1.48 1.38 1.28 1.18 1.08 0.98 0.88 0.78 0.68 0.58 xy -0.5 0 0.5-0.200.20.40.6  0.94 0.84 0.74 0.64 0.54 0.44 0.34 0.24 0.14 34 (a) (b) Fig. 17 a) pressure coefficient distribution over the hydrofoil surfaces and b) their corresponding shape and the cavity pocket indicated by 𝐶𝑝=−𝜎. x-Cp -0.4 -0.2 0 0.2 0.4-1-0.500.51 Optimized-fronthydrofoil Optimized-rearhydrofoil Original x-0.6 -0.4 -0.2 0 0.2 0.4 0.6 35 6. Conclusion In this paper , a slotted hydrofoil is optimized to find out the slot location and the slot angle that results in maximizing the lift coefficient and minimizing the total length of the cavity . Through the optimization procedure, a computational fluid dynamic (CFD) solver is required to provide the lift coefficient and the total length of the cavity for the optimizer. For this aim, the Frink numerical method is newly developed here to discretize the spatial terms in the preconditioned Euler equations based on the artificial compressibility method. The artificial dissipation terms are added to the se equations for stabili ty reasons and a sensor is used to retain the solution accuracy in smooth regions . Moreov er, source terms that model the cavitation process are also added to the equations. The optimizer is constructed using a surrogated model based on the deep ANN surrogate model and the genetic algorithm. The Latin Hypercube algorithm is used to generate initial sample points in the design space . The outcome of the optimization process is an optimized slotted NACA 66 (MOD) hydrofoil satisfying the cost functions defined. Some conclusions and remarks regarding the present work are as follows: - The comparison between results obtained by the CFD solv er and available experimental results shows that the Frink finite volume -based solver newly proposed and developed here can be effectively used for simulating cavitating flows and can be considered as an alternative for the other cavitating flow solvers. - Regarding the original hydrofoil NACA 66 (MOD) , one cavity pocket is formed on its suction side. However, in the slotted hydrofoil, two cavities will emerge on the hydrofoil surface, namely, the leading -edge and mid -chord cavities. Smaller leading -edge cavity size will be found around 𝜃=45° and at high slot angles. In the case of the mid-chord cavity, higher slot angles should be used with a slot placed around the middle . - From the results obtained it can be concluded that high slot angles very close to the leading edge desirably result in the small total length of the cavity 𝐿𝑐 and high lift coefficient 𝐶𝑙. - The deep ANN surrogated mode is examined against CFD results wh ich shows a good agreement with an error lower than 4% , thus it can be reliably used in the genetic algorithm for optimizing purposes. - Results obtained for the optimized slotted hydrofoil show that with a penalty of about 11.4% in the lift coefficient, it is possible to decrease the cavity size by about 58.4%. 36 - The contours for the lift coefficient , the total length of the cavity, the length of the leading - edge cavity, and the length of the mid -chord cavity obtained in this study can be used by engineers to figure out how deviation from the optimum point can affect the performance of their artifacts . Acknowledgment Mahya Hajihassanpour would like to thank Prof. Kazem Hejranfar and Dr. Amir Hossein Khodabakhsh for their technical advice .