An Adjoint Optimization Prediction Method for Partially Cavitating Hydrofoils

: Much work has been done over the past years to obtain a better understanding, predict and alleviate the effects of cavitation on the performance of lifting surfaces for hydrokinetic turbines and marine propellers. Lifting-surface sheet cavitation, when addressed as a free-streamline problem, can be predicted up to a desirable degree of accuracy using numerical methods under the assumptions of ideal ﬂow. Typically, a potential solver is used in conjunction with geometric criteria to determine the cavity shape, while an iterative scheme ensures that all boundary conditions are satisﬁed. In this work, we propose a new prediction model for the case of partially cavitating hydrofoils in a steady ﬂow that treats the free-streamline problem as an inverse problem. The objective function is based on the assumption that on the cavity boundary, the pressure remains constant and is evaluated at each optimization cycle using a source-vorticity BEM solver. The attached cavity is parametrized using B-splines, and the control points are included in the design variables along with the cavitation number. The sensitivities required for the gradient-based optimization are derived using the continuous adjoint method. The proposed numerical scheme is compared against other methods for the NACA 16-series hydrofoils and is found to predict well both the cavity shape and cavitation number for a given cavity length.


Introduction
One of the most demanding technical requirements imposed on the design of lifting surfaces for nozzles, water pumps, hydrokinetic turbines, and marine propulsion systems is due to cavitation [1,2]. As it develops, cavitation creates noise, vibration, metal erosion, and a drop in performance [3,4]. Over the past years, researchers have obtained a better understanding of cavitation due to the standardization of experimental research, intending to provide an accurate experimental database for the validation of computational methods [5,6]. Moreover, numerical predictions can provide useful information and details about the onset and evolution of cavitation. Therefore, simulating various forms of cavitation using multiphase CFD has become increasingly common in the last several years [7][8][9][10].
CFD-based design optimization is a relatively new and emerging field that provides a direct link between numerical simulations and the required design improvements. In the work of [11], a continuous adjoint method is developed for the design optimization of a cavitating hydrofoil based on a homogeneous multiphase mixture model. The continuous adjoint method is also employed in the work of [12] for the shape optimization of hydraulic turbomachines, with objective functions targeting various aspects of design improvements. However, applications of the continuous adjoint method for cavitation problems are not limited to CFD solvers. This method based on an ideal flow solver was successfully employed by the authors of [13] to determine the optimum shape of a supercavitating torpedo in terms of drag minimization given certain operating conditions, along with through a convergence study and comparisons against finite differences (FDM). The effects of hydrofoil thickness on the cavitation number and cavity volume are investigated in a parametric case study for hydrofoils based on the NACA 16-series. Finally, concluding remarks are provided along with suggestions for future work and research directions.

Inverse Problem
This section presents the mathematical formulation behind the proposed adjoint optimization prediction method for the case of partially cavitating hydrofoils in a steady flow. The formulation of the inverse problem (see Figure 1) and the assumptions made are discussed in detail. Particularly, the governing partial differential equations for the lifting flow problem of a hydrofoil in steady flow with a given sectional profile and inflow conditions, namely the primal problem, can be solved numerically in the sense of boundary integral equations (BIE), as discussed in Section 2. 1. tigated in a parametric case study for hydrofoils based on the NACA 16-series. Finally, concluding remarks are provided along with suggestions for future work and research directions.

Inverse Problem
This section presents the mathematical formulation behind the proposed adjoint optimization prediction method for the case of partially cavitating hydrofoils in a steady flow. The formulation of the inverse problem (see Figure 1) and the assumptions made are discussed in detail. Particularly, the governing partial differential equations for the lifting flow problem of a hydrofoil in steady flow with a given sectional profile and inflow conditions, namely the primal problem, can be solved numerically in the sense of boundary integral equations (BIE), as discussed in Section 2.1.
In the present work, however, we address the inverse problem of determining the cavity shape and cavitation number for a given cavity length and inflow conditions. The optimal solution must coincide with a streamline, and therefore, we introduce the primal problem as a constraint for the inverse problem to ensure that the optimal solution is found within the context of admissible solutions. The sensitivities required for the first order gradient-based optimization scheme are obtained using the continuous adjoint method. The adjoint boundary value problem (BVP) is also solved numerically in the sense of BIE using the same numerical method as the primal. Regarding the problem statement, the "fixed cavity length" simplification facilitates the verification of the proposed numerical scheme through comparisons against other methods found in literature, i.e., the work of [19,22], that follow the same assumptions. However, the "fixed cavity length" assumption can be waived by means of an iterative scheme, as shown in the work of [19], in order for our method to consider more realistic applications in future work. For the parametric representation of the hydrofoil with the attached cavity (see Figure 1), we adopt a clock-wise convention, In the present work, however, we address the inverse problem of determining the cavity shape and cavitation number for a given cavity length and inflow conditions. The optimal solution must coincide with a streamline, and therefore, we introduce the primal problem as a constraint for the inverse problem to ensure that the optimal solution is found within the context of admissible solutions. The sensitivities required for the first order gradient-based optimization scheme are obtained using the continuous adjoint method. The adjoint boundary value problem (BVP) is also solved numerically in the sense of BIE using the same numerical method as the primal.
Regarding the problem statement, the "fixed cavity length" simplification facilitates the verification of the proposed numerical scheme through comparisons against other methods found in literature, i.e., the work of [19,22], that follow the same assumptions. However, the "fixed cavity length" assumption can be waived by means of an iterative scheme, as shown in the work of [19], in order for our method to consider more realistic applications in future work. For the parametric representation of the hydrofoil with the attached cavity (see Figure 1), we adopt a clock-wise convention, with x(t o ) = x(t 1 ), y(t o ) = y(t 1 ) denoting the trailing edge (TE). The hydrofoil/cavity outline (highlighted in gray and blue in Figure 1) can be re-modeled for any given set of nodal coordinates using interpolation, with more details presented in Section 3.1 that follows.
To determine the initially unknown shape of the cavity for a given cavity length, additional information regarding the cavity termination region is essential to the modeling. In this work, we implement the cavity termination model presented in the work of Kinnas et al. [19]. Particularly, the detachment and termination points of the constantpressure cavity region are denoted by s D , s T respectively, whereas s L is the reattachment point as shown in Figure 1. Based on the parametric representation of the curve r dt with {t o , t D , t T } ∈ I. The dot notation . r(t) denotes differentiation with respect to the curve parameter t.
In the context of the primary equations, the attached cavity shape is based on an initial guess. The methodology for determining the optimal shape of the cavity boundary Γ c = (x, y) ∈ R 2 : t ∈ (t D , t T ) that is presented here refers to an optimization process that produces improved cavity shape estimates (i.e., the boundary segment highlighted in blue in Figure 1) at each optimization cycle. On the unknown boundary, we assume that the cavitation number is constant and the pressure uniform, an assumption that yields the following cost functional, where p υ is the previously unknown target (vaporization) cavity pressure and p the pressure estimate obtained from the primal solver, relying solely on the outline of the hydrofoil/cavity and the inflow conditions. The cavitation number is defined as where U ∞ , p ∞ and T ∞ are, respectively, the reference velocity, pressure, and temperature in the flow (usually upstream quantities), ρ the density of the fluid, and p υ (T ∞ ) is the saturated vapor pressure. Particularly, the kernel of the objective function can also be expressed in terms of the non-dimensional pressure coefficient (C p = 0.5ρU 2 ∞ ) as follows Notably, admissible solutions in terms of design variables b n , n = {1, . . . , N} must comply with the requirements of incompressible, inviscid, and irrotational fluid motion in the region, i.e., the work of [24,25]. The design variable vector consists of the unknown cavitation number and parameters that affect the shape of the attached cavity, to be discussed in Section 3.1. The primal BVP; see also Section 2.1, with respect to the disturbance velocity potential ϕ serves as a constraint for the optimization problem, formulated as with the Kutta condition in Equation (6) expressing the assumption that the tangential velocity on the upper and lower sides of the TE when superimposed is equal to zero; thus, the velocity in the vicinity of the TE is finite. It is important to note that the primal variable ϕ = ϕ(x, y; b n ) is implicitly dependent on the design variables. The above equations are paired with the following "no-entry" (zero flux) boundary condition wheren,τ denote the unit normal and tangential vectors on the boundary of the hydrofoil. Finally, all admissible solutions must also comply with a condition at infinity Note that in accordance with other works found in the literature, such as the work of [19], on the transition region (x, y) ∈ R 2 : t ∈ (t T , t L ) only the "no-entry" boundary condition given in Equation (7) is to be satisfied. A straightforward approach to obtain an estimate of the sensitivity derivatives for each design variable δF/δb n would be to implement finite differences. However, for central finite differences, this approach requires 2N evaluations of the primal solver, with N denoting the total number of design variables. As an alternative, we implement the continuous adjoint method to produce estimates of the sensitivity derivatives that require fewer evaluations of the primal solver and only two evaluations per optimization cycle.
The derivation of sensitivities for the continuous adjoint method occurs analytically at the level of the partial differential equations. A standard approach for the derivation of the adjoint-state equations is via a Lagrange multiplier ψ, denoted as the adjoint velocity potential, that is continuous and twice differentiable; see the work of [11,12,23]. It is worth mentioning here that the present formulation is based on the BVP in Equations (5)- (8) and the minimization of the cost functional F as defined by Equation (2), whereas the numerical solution of the primal and adjoint BVPs is obtained in the sense of BIE as discussed in Section 2.1 that follows. Instead of using analytical techniques, the derivation of the sensitivities can be based on symbolic mathematics to further facilitate the process, as shown in the work of [31], where an adjoint solver based on the classical Douglas-Neumann panel method is used to tackle the design optimization problem of airfoil sections.
The augmented cost functional consists of integrals on the flow domain Ω and its boundary consisting of the original hydrofoil geometry superimposed with the attached cavity ∂Ω = Γ c ∪ Γ w , where δ(s − s o ) is the Dirac delta generalized function. To obtain the sensitivity derivatives, we take the first variation of the functional, i.e., the work of [32], with respect to the design variables b n , n = {1, . . . , N} as follows The first term in Equation (10) becomes The first two terms in Equation (11) appear after the implementation of the Leibnitz rule of integration. The differentiated form of the cost functional is, therefore, Note that it is not possible to determine δϕ/δb n explicitly since the velocity potential is implicitly dependent on the design variables. Using Gauss theorem in conjunction with the boundary condition in Equation (7), the second term in Equation (9) yields The third term in Equation (9), imposing a stagnation point at the vicinity of the TE, becomes The tangential vectors in the vicinity of the TE are not dependent on the design variables and therefore δτ o /δb n = δτ 1 /δb n = 0. The first variation of the augmented functional is derived from Equations (10)- (14), after re-arranging terms and implementing integration by parts, where serves as the approximation formula for the sensitivity derivatives. The adjoint boundary value problem is introduced so that Equation (16) becomes independent of δϕ/δb n . The boundary value problem is very similar to the primal; however, here on the cavity boundary, a non-homogeneous Neumann condition is to be satisfied Since the admissible cavity shapes have fixed end-points s D , s L , the variation at these positions along the boundary is zero δϕ/δb n . In the present work, the pressure attenuation region t ∈ (t D , t T ) is defined by numerical testing, and it is found to be compatible with similar findings in the literature, i.e., the work of [18,19]. Based on our formulation at the beginning of the pressure attenuation region, the velocity must be equal to V t = U ∞ √ σ + 1, or equivalently equal to p = p υ . The undefined constant associated with the exterior Neumann boundary value problem Equations (17)- (19) is resolved by assuming that ψ = 0 at the TE. The latter is used in conjunction with the continuity assumption of the adjoint potential, thus allowing us to drop the last two terms on the right-hand side of Equation (15).

Representation of the Solution Field
The potential ϕ and velocity field V = ∇ϕ, corresponding to the fluid flow problem, are the solutions to the BVP presented in Equations (5)- (8). This solution must also satisfy the equivalent weak form of the BVP in the sense of BIE. One standard approach to derive the integral equation is via Green's theorem based on the fundamental solution of the governing equation and the perturbation potential ϕ as the targeted unknown, i.e., the work of [24,25]. However, in the present work, we follow an approach based on representing the unknown perturbation velocity field by superimposing solutions of the governing equation where G 1 (r, r o ) and G 2 (r, r o ) denote the point source-and vortex-type singularities, solutions to the Laplace equation that also satisfy by definition the required condition at infinity in Equation (8). We differentiate Equation (20) under the integral sign and use a limiting process to derive the integral equation that holds for each point on the boundary [25]. Then the no-entry boundary condition is formulated in the sense of BIE as follows Working with the BIE yields a significant dimensionality reduction since it is possible to determine the unknown velocity field by solving Equation (23), which is solely dependent on boundary data. The targeted unknowns for the integral equation are the strengths of the source and vortex-type singularities σ(s), γ.

Geometry Parametrization
For the NACA 16-airfoil series, the symmetric thickness distribution is defined using two concatenated cubic Bezier curves as x f oil , y f oil , based on the work of [33], thus allowing the implementation of cosine spacing near the leading/trailing edge. To obtain the hydrofoil outline with the attached cavity that leads to the parametric representation in Equation (1), we used a parabolic function with one free parameter β ∈ R + that controls the initial cavity volume to generate an initial guess. Based on the "fixed cavity length" assumption, the initial guess for the attached cavity can be defined on the portion of the boundary where t ∈ (t D , t L ). On this part of the boundary, we superimpose the original hydrofoil section with the following.
The set of points generated from Equation (24) is then used to define the parametric representation of the attached cavity based on B-spline interpolation of the fourth order.
In this work, we chose B-splines for their characteristic locality property [27] since a perturbation of the control points of a Bezier curve affects the entire shape of the curve. The B-spline control points that determine the attached cavity geometry define the set of design variables b n , n = {1, . . . , N} along with the unknown cavitation number. The initial guess significantly affects the convergence rate of the optimization method, and this aspect is discussed in detail in Section 4. Embedding such detail into the geometric model that generates the initial guess for the attached cavity is quite useful when simulating cavitating hydrofoils at various cavity lengths, as it affects the convergence rates. Starting from a parabolic shape does not guarantee smooth transition near the cavity closure region, and therefore we must also adjust the slope of the B-spline at the vicinity of the cavity termination region, as shown in Figure 2. The same holds for the attached cavity in the vicinity of the leading edge (LE) to ensure that the original hydrofoil geometry is not compromised by the cavity boundary during the optimization process.
region. As previously discussed, the extent of the pressure attenuation region TL is an input parameter for the proposed model, as adopted by the authors of [19]. An interesting discussion about the pressure attenuation region and its effects on the pressure profile is found in the work of [34]. It is also important to note that, based on this geometric representation, the wall termination model (see the work of [22]) can be easily implemented as an alternative model for the cavity termination region. Since our optimization is subjected to no constraints, it is essential to verify during the optimization process that the attached cavity remains compatible with the physical assumption of the free-streamline problem, i.e., the attached cavity boundary does not intersect with the original hydrofoil boundary. The only intersections allowed are the detachment and reattachment points of the cavity. The first and last two control points (see Figure 2) are not included in the design variable vector.

Solution of the Primal Problem
For the solution of the primal problem in Equation (23), a low order panel method based on piece-wise constant source and vortex distributions is implemented; see the work of [24,25]. The boundary is approximated with 1 M straight line segments, known as panel elements For the optimization problem, the design variables consist of the cavitation number, the coordinates of the B-spline control points, and the length of the pressure attenuation region. As previously discussed, the extent of the pressure attenuation region λ TL is an input parameter for the proposed model, as adopted by the authors of [19]. An interesting discussion about the pressure attenuation region and its effects on the pressure profile is found in the work of [34]. It is also important to note that, based on this geometric representation, the wall termination model (see the work of [22]) can be easily implemented as an alternative model for the cavity termination region.
Since our optimization is subjected to no constraints, it is essential to verify during the optimization process that the attached cavity remains compatible with the physical assumption of the free-streamline problem, i.e., the attached cavity boundary does not intersect with the original hydrofoil boundary. The only intersections allowed are the detachment and reattachment points of the cavity. The first and last two control points (see Figure 2) are not included in the design variable vector.

Solution of the Primal Problem
For the solution of the primal problem in Equation (23), a low order panel method based on piece-wise constant source and vortex distributions is implemented; see the work of [24,25]. The boundary is approximated with M 1 straight line segments, known as panel elements The total velocity V i = (u i , υ i ) at the midpoint of the i-th panel is given by where r si = (x si , y si ) is the relative position of the control point (x i , y i ) to the point of integration, particularly x si = x i − x s , y si = y i − y s . The discretized form of the boundary conditions is as follows Based on the above, Equations (28) and (29) form a linear system of M 1 + 1 equations with M 1 + 1 unknowns with respect to the strengths of the source and vortex distributions. The pressure coefficient at the midpoint of each panel can be calculated using After a solution has been obtained, we can verify that the results in terms of lift CL and moment coefficient CM are in suitable agreement with other methods found in the literature. In Figure 3a, we present a comparison between the present panel method (BEM), experimental data from the work of [35], and inviscid simulations performed with XFOIL [36] for the case of a NACA 4412 airfoil at U ∞ = (1, 0) m/s inflow. The present method is found to be in excellent agreement with the inviscid simulations and to compare well against experimental data, considering the limitations of an ideal flow-based method. In Figure 3b, we showcase that the primal solver converges to the solution of the finest discretization. This test case reveals that a spatial discretization of 400 panels gives lift and moment coefficient predictions below the threshold of 1% error; thus, this spatial discretization is selected for the test cases presented in the following sections.
Based on the above, Equations (28) and (29)  1 , After a solution has been obtained, we can verify that the results in terms of lift CL and moment coefficient CM are in suitable agreement with other methods found in the literature. In Figure 3a, we present a comparison between the present panel method (BEM), experimental data from the work of [35], and inviscid simulations performed with XFOIL [36] for the case of a NACA 4412 airfoil at (1,0) m/s U inflow. The present method is found to be in excellent agreement with the inviscid simulations and to compare well against experimental data, considering the limitations of an ideal flow-based method. In Figure 3b, we showcase that the primal solver converges to the solution of the finest discretization. This test case reveals that a spatial discretization of 400 panels gives lift and moment coefficient predictions below the threshold of 1% error; thus, this spatial discretization is selected for the test cases presented in the following sections.

Solution of the Adjoint Problem
The adjoint field equation accompanied with a Neumann-type boundary condition can also be treated in the sense of BIE. However, the numerical solution of this problem is much simpler than the primal problem. The same discretization scheme is used to approximate the boundary of the foil/cavity. Since this problem does not have circulation, vortex elements are not used in the representation, and on each panel, we simply have a piece-wise constant distribution of sources, σ(s) = σ j , j = {1, . . . M 1 }. The adjoint velocity V ai = (u ai , υ ai ) at the midpoint of the i-th panel, taking into account the method of images, is given by with u sij , υ sij defined in Equation (27). The boundary conditions in the discretized form are In a similar manner, Equation (31) forms a linear system of equations with respect to the unknown strengths of source distributions σ j .
To illustrate the qualitative differences between the primal and the adjoint solutions, we present in Figure 4 results obtained with the present method, for the case of a NACA16-006 hydrofoil at a = 4 • , U ∞ = (1, 0) m/s. The cavity shape in this example is an initial guess for the test case to be presented in the section that follows. Particularly, the cavity length is set equal to l c /c = 0.5. In Figure 4a, we present in color-map format the non-dimensional pressure field superimposed with the disturbance velocity streamlines as obtained from the fully wetted algorithm (primal problem). In Figure 4b, we present the corresponding adjoint potential field superimposed with the adjoint velocity streamlines. ages, is given by 11 , , NN ai j sij ai j sij jj uu (31) with ,  (32) In a similar manner, Equation (31) forms a linear system of equations with respect to the unknown strengths of source distributions j .
To illustrate the qualitative differences between the primal and the adjoint solutions, we present in Figure 4 results obtained with the present method, for the case of a NACA16-006 hydrofoil at 4deg a , (1,0) m/s U . The cavity shape in this example is an initial guess for the test case to be presented in the section that follows. Particularly, the cavity length is set equal to / 0.5 c lc . In Figure 4a, we present in color-map format the non-dimensional pressure field superimposed with the disturbance velocity streamlines as obtained from the fully wetted algorithm (primal problem). In Figure 4b, we present the corresponding adjoint potential field superimposed with the adjoint velocity streamlines.

Calculation of Sensitivity Derivatives
The integrals in Equation (16) are evaluated by means of standard numerical integration, and the derivatives of the geometric quantities, such as the normal/tangential vectorsn,τ or the Jacobian . r , are approximated with a central FDM scheme based on information obtained from the B-spline representation. The sensitivities obtained with the present method, denoted as δL/δb n , are compared against results obtained with the direct implementation of central differences on the cost functional δF/δb n for verification. It is important to note that calculating the sensitivities with central finite differences requires 2N evaluations of the BEM solver, with N denoting the total number of design variables. On the other hand, the method we propose for approximating the sensitivities based on the continuous adjoint method requires only one additional evaluation of the BEM solver.

The Minimization Algorithm
The simple steepest-descend method is implemented to treat numerically the deterministic optimization problem. Our analysis leads to the PCavPreMod algorithm described below in Algorithm 1. The rate of convergence is strongly dependent on the selection of the step size η, as well as the initial guess for the design variable vector. It is important to note that since the step size η is dimensional, the convergence thresholds must be properly tuned beforehand. Assuming that the constants selected for the steepest-descend scheme are suf-ficiently small, a local optimum is reached as the design cycles progress. The optimization process continues until the solution has converged, i.e., the values of the sensitivities and the cost functional are below a certain threshold, or the number of maximum evaluations M max has been reached. One way to explore further the design space within the capabilities of a gradient-based method is to search for the optimal solution starting from different initial cavity shapes. Solve primal problem to define ϕ Solve adjoint problem to define ψ Calculate sensitivitie derivatives SD = δL/δb n Update b k+1

Verification of the Present Method
In the context of model verification, two study cases were selected from the work of [18], where the numerical simulation of the flow around two-dimensional partially cavitating hydrofoils is also based on the assumptions of ideal flow. For both cases, a constant cavity length is considered (l c /c = 0.5), whereas the cavitation number and cavity shape are to be determined upon the solution of the inverse problem presented in Section 2. These study cases refer to a NACA16-006 section at a = 4 • and a = 6 • angle of attack. The pressure attenuation region is taken as λl c = 0.1l c for both cases. The initial guess for the cavity boundary for all simulations is based on Equation (24) and s D is positioned at the LE. The first step toward the verification of the present method is to achieve suitable agreement between the sensitivities obtained with the adjoint method and finite differences. It is also important to show that convergence to an optimal solution is achieved as the optimization cycles progress.
In Figure 5a, we present the sensitivities for each design variable or degree of freedom (DOF), i.e., the cavitation number, the yi-and xi-coordinates of the B-spline control points. The optimization algorithm requires information about the step size, a vector equal in size with the design variables vector. It is possible for the designer to restrict the xi-coordinates of the B-spline cavity by simply setting the corresponding step sizes equal to zero. As we discussed in a previous section about the geometric representation of the cavity, the optimization algorithm we implemented is unconstrained, and therefore, it is important that the step sizes that correspond to the xi-coordinates, especially in the vicinity of the LE, are set to zero to ensure that the geometry of the hydrofoil is not compromised. On the other hand, the control points near the cavity termination region must be allowed to have the freedom to move with respect to both coordinate directions, so the step sizes of these design variables are tuned accordingly. In Figure 5b, we observe the evolution of the cost functional values as the optimization cycles progress for NACA16-006 at a = 4 • , with a clear indication that there exists a number of cycles where the cost functional convergences to values lower than 1 × 10 −4 .
The convergence of the present method is dependent on the cavity initialization since it affects the sensitivity derivatives, especially near the edges of the cavity. In order to illustrate this effect in Figure 6a, we present three different initializations for the attached cavity, namely the "h1", "h2", and "h3" shapes. These initial guesses correspond to values β = {1.3, 1.2, 0.5} of the parameter in Equation (24). In Figure 6b, we present the evolution of the cost functional with respect to the optimization cycle for each initialization. All other parameters are kept the same. It is seen that all three initializations converge, and the results concerning the cavitation number σ and cavity volume prediction Vol/c 2 are provided in Table 1. Suitable selections for the starting shape of the cavity can be determined by the systematic application of the present method and numerical experimentation with various foil geometries and inflow conditions, and this is left for future work. To further highlight the ability of the proposed numerical scheme to predict both the cavity shape and the cavitation number as part of the solution, we present in Figures 7 and 8 the pressure coefficient distribution as obtained with the present method in comparison with results from the work of [18].
the optimization cycles progress.
In Figure 5a, we present the sensitivities for each design variable or degree of freedom (DOF), i.e., the cavitation number, the yi-and xi-coordinates of the B-spline control points. The optimization algorithm requires information about the step size, a vector equal in size with the design variables vector. It is possible for the designer to restrict the xicoordinates of the B-spline cavity by simply setting the corresponding step sizes equal to zero. As we discussed in a previous section about the geometric representation of the cavity, the optimization algorithm we implemented is unconstrained, and therefore, it is important that the step sizes that correspond to the xi-coordinates, especially in the vicinity of the LE, are set to zero to ensure that the geometry of the hydrofoil is not compromised. On the other hand, the control points near the cavity termination region must be allowed to have the freedom to move with respect to both coordinate directions, so the step sizes of these design variables are tuned accordingly. In Figure 5b, we observe the evolution of the cost functional values as the optimization cycles progress for NACA16-006 at a = 4deg, with a clear indication that there exists a number of cycles where the cost functional convergences to values lower than The convergence of the present method is dependent on the cavity initialization since it affects the sensitivity derivatives, especially near the edges of the cavity. In order to illustrate this effect in Figure 6a, we present three different initializations for the attached cavity, namely the "h1", "h2", and "h3" shapes. These initial guesses correspond to values {1.3, 1.2, 0.5} of the parameter in Equation (24). In Figure 6b, we present the evolution of the cost functional with respect to the optimization cycle for each initialization. All other parameters are kept the same. It is seen that all three initializations converge, and the results concerning the cavitation number and cavity volume prediction 2 / Vol c are provided in Table 1. Suitable selections for the starting shape of the cavity can be determined by the systematic application of the present method and numerical experimentation with various foil geometries and inflow conditions, and this is left for future work.
To further highlight the ability of the proposed numerical scheme to predict both the cavity shape and the cavitation number as part of the solution, we present in Figures 7  and 8 the pressure coefficient distribution as obtained with the present method in comparison with results from the work of [18].  The initial and final cavity shapes are provided as well. The results presented in Figure 8 correspond to the "h2" initial cavity that showed the best convergence. These results were directly compared with the method originally proposed in the work of [19] but also with CFD simulations for partially cavitating hydrofoils in a steady flow. The pressure profiles obtained with potential-based methods, such as ours, used for the prediction of sheet cavitation are in accordance with the CFD simulations apart from the transition zone region, where the closure model is used [19]. In Figures 7a and 8a, we can also observe that for higher angles of attack, the cavity sectional area increases.   The initial and final cavity shapes are provided as well. The results presented in Figure 8 correspond to the "h2" initial cavity that showed the best convergence. These results were directly compared with the method originally proposed in the work of [19] but also with CFD simulations for partially cavitating hydrofoils in a steady flow. The pressure profiles region, where the closure model is used [19]. In Figures 7a and 8a, we can also observe that for higher angles of attack, the cavity sectional area increases.
(a) (b) Figure 8. Comparison between the present method for the results found in the work of [18]. (a) Initial (h2) and final cavity shape for NACA16-006 with lc = 0.5 at a = 6deg. The squares represent the control points of the B-spline representation. (b) Pressure coefficient.

Effects of Thickness on Cavity Volume and Cavitation Number
The cavity volume, or more precisely, the sectional area of the cavity, is also of great interest since, in an unsteady case, changes in volume determine the monopole-type acoustic source strength; see [3,4]. To investigate the effects of thickness on the sectional area of the cavity, we performed simulations for the NACA16-006, 16-009, and 16-012 sections at a = 4deg angle of attack. This allows a direct comparison between the present method and the results obtained with the surface singularity method introduced by Uhlman in the work of [22]. The predicted sectional area of the cavity is presented in Fig-Figure 8. Comparison between the present method for the results found in the work of [18]. (a) Initial (h2) and final cavity shape for NACA16-006 with lc = 0.5 at a = 6 • . The squares represent the control points of the B-spline representation. (b) Pressure coefficient.

Effects of Thickness on Cavity Volume and Cavitation Number
The cavity volume, or more precisely, the sectional area of the cavity, is also of great interest since, in an unsteady case, changes in volume determine the monopole-type acoustic source strength; see [3,4]. To investigate the effects of thickness on the sectional area of the cavity, we performed simulations for the NACA16-006, 16-009, and 16-012 sections at a = 4 • angle of attack. This allows a direct comparison between the present method and the results obtained with the surface singularity method introduced by Uhlman in the work of [22]. The predicted sectional area of the cavity is presented in Figure 9 as a function of the cavity length. The numerical parameters that correspond to the simulations performed for NACA16-006 hydrofoil with the adjoint optimization prediction method are given in Table 2 below, with l c denoting the cavity length, N DOFS the total number of DOFs, λ% the extent of the transition region as a percentage of the cavity length, and Ncycles the maximum number of evaluations to achieve convergence. The values of parameter β used to generate each initial guess are also included in Table 2. It is important to note that similar values are used for the other hydrofoil sections. The pressure profiles and predicted cavity for a NACA16-006 hydrofoil at a = 4 • for cavity lengths l c = {0.2, 0.4, 0.6, 0.8} are shown in Figure 10.
To verify that the proposed method is capable of predicting the sensitivity derivatives, we compared results obtained with the present method against central differences for a selected study case from the work of [18], as shown in Figure 5a. The adjoint sensitivities are found to be in suitable agreement with the FDM both in terms of magnitude and sign. These results also justify the additional assumptions made in the mathematical modeling in Section 2. The PCavPreMod algorithm for both cases converges to an optimal solution, as indicated by the results presented in Figures 5b and 6b. The cavity initialization affects not only the sensitivities near the edges of the cavity but also the convergence rate of the method as shown in Figure 6b. It is also important to note that the steepestdescend method has a relatively slow convergence among the various first order gradientbased optimization algorithms.
The benefits of using the continuous adjoint method to compute the sensitivities are illustrated in these comparisons, since computing the sensitivities with the present method requires only two evaluations of the BEM solver, i.e., one for the primal problem and one for the adjoint, regardless of the number of design variables. For the same case, approximating the sensitivities using central differences would require two evaluations per design variable and a total of 2N evaluations. The proposed method makes it also possible for the designer to restrict some of the design variables during the early optimization cycles, for example, the xi-coordinates of the control points near the LE, and to restore them by changing the step size later on to achieve a better pressure profile without compromising the computational cost of each optimization cycle. Gradient-based optimization with FDM is not that flexible and most certainly computationally intensive. Figure 9. Comparison between the present method (PM) and numerical results found in Uhlman [22] for the NACA 16-006, 16-009, 16-012 sections at a = 4deg. The sectional area is presented as a function of cavity length. Figure 9. Comparison between the present method (PM) and numerical results found in Uhlman [22] for the NACA 16-006, 16-009, 16-012 sections at a = 4 • . The sectional area is presented as a function of cavity length.  In Figures 7 and 8, we proceed by presenting the pressure profiles and corresponding cavity shapes for the NACA16-006 at a = 4deg and a = 6deg, respectively. The pressure profiles and cavitation number predictions are found to be in accordance with the numerical results published in the work of [18]. However, it is important to note that a wave-like behavior of the pressure profile is observed near the cavity detachment and reattachment points. This does not seem to affect the prediction of the cavitation number or the cavity sectional area, as shown in Figure 9. One explanation for this behavior comes from the definition of the cost functional itself in Equation (2), which does not contain any information about the derivative of the pressure profile. This issue could be resolved by adding terms that contain the derivatives of the pressure coefficient in the cost functional.
So far, the proposed model has been verified in terms of predicting the pressure profile and cavitation number for cavity lengths equal to half the chord. The parametric case study in Figure 9 reveals that the developed numerical scheme is also suitable for the prediction of the cavity's sectional area for a wide range of cavity lengths. This case study is taken from the work of Uhlman [22]. A specific set of parameters needs to be tuned prior to each simulation contained in Figure 9. Some of these parameters are included in Table  2 for the case of a NACA16-006 at a = 4deg angle of attack. Since our method is based on an optimization principle, similar results can be obtained with minor changes to these parameters; nevertheless, some interesting trends are observed. The length of the pressure attenuation region, or transition region, is between 0.05% and 0.22% as a percentage of the

Discussion
To verify that the proposed method is capable of predicting the sensitivity derivatives, we compared results obtained with the present method against central differences for a selected study case from the work of [18], as shown in Figure 5a. The adjoint sensitivities are found to be in suitable agreement with the FDM both in terms of magnitude and sign. These results also justify the additional assumptions made in the mathematical modeling in Section 2. The PCavPreMod algorithm for both cases converges to an optimal solution, as indicated by the results presented in Figures 5b and 6b. The cavity initialization affects not only the sensitivities near the edges of the cavity but also the convergence rate of the method as shown in Figure 6b. It is also important to note that the steepest-descend method has a relatively slow convergence among the various first order gradient-based optimization algorithms.
The benefits of using the continuous adjoint method to compute the sensitivities are illustrated in these comparisons, since computing the sensitivities with the present method requires only two evaluations of the BEM solver, i.e., one for the primal problem and one for the adjoint, regardless of the number of design variables. For the same case, approximating the sensitivities using central differences would require two evaluations per design variable and a total of 2N evaluations. The proposed method makes it also possible for the designer to restrict some of the design variables during the early optimization cycles, for example, the xi-coordinates of the control points near the LE, and to restore them by changing the step size later on to achieve a better pressure profile without compromising the computational cost of each optimization cycle. Gradient-based optimization with FDM is not that flexible and most certainly computationally intensive.
In Figures 7 and 8, we proceed by presenting the pressure profiles and corresponding cavity shapes for the NACA16-006 at a = 4 • and a = 6 • , respectively. The pressure profiles and cavitation number predictions are found to be in accordance with the numerical results published in the work of [18]. However, it is important to note that a wave-like behavior of the pressure profile is observed near the cavity detachment and reattachment points. This does not seem to affect the prediction of the cavitation number or the cavity sectional area, as shown in Figure 9. One explanation for this behavior comes from the definition of the cost functional itself in Equation (2), which does not contain any information about the derivative of the pressure profile. This issue could be resolved by adding terms that contain the derivatives of the pressure coefficient in the cost functional.
So far, the proposed model has been verified in terms of predicting the pressure profile and cavitation number for cavity lengths equal to half the chord. The parametric case study in Figure 9 reveals that the developed numerical scheme is also suitable for the prediction of the cavity's sectional area for a wide range of cavity lengths. This case study is taken from the work of Uhlman [22]. A specific set of parameters needs to be tuned prior to each simulation contained in Figure 9. Some of these parameters are included in Table 2 for the case of a NACA16-006 at a = 4 • angle of attack. Since our method is based on an optimization principle, similar results can be obtained with minor changes to these parameters; nevertheless, some interesting trends are observed. The length of the pressure attenuation region, or transition region, is between 0.05% and 0.22% as a percentage of the cavity length. The number of design variables changes between 11 and 30, with the smaller values corresponding to the smallest cavity lengths. This is reasonable since overfitting a B-spline curve is more prone to these wave-line pressure profiles observed in Figures 7b and 8b. The number of optimization cycles required to achieve convergence after proper selection of the steepest-descend step size was between 250 and 300 iterations.
In Figure 10, we present the pressure profiles corresponding to the NACA16-006 at a = 4 • for various cavity lengths l c = {0.2, 0.4, 0.6, 0.8}. The pressure profile fluctuations are mostly concentrated in the LE region, whereas the pressure profile at the transition zone is very similar to the cavity termination model proposed by Kinnas et al. [18], as was expected.

Conclusions
The problem of steady, partially cavitating two-dimensional hydrofoils with known cavity length is addressed as an inverse problem. The continuous adjoint method based on an ideal flow model is used for the derivation of the sensitivities required for first-order gradient-based optimization. The attached cavity is parametrized using B-splines, and the control points are included in the set of design variables along with the unknown cavitation number. The proposed numerical scheme is compared with other methods that follow the "fixed cavity length" assumption and is found to predict well the cavity shape, volume, and cavitation number. The benefits of using the adjoint method to predict the sensitivity derivatives instead of FDM are illustrated in a selected study case. The proposed method is also found to predict well the effects of thickness on the sectional area of the cavity based on a parametric case study for the NACA 16-006, 16-009, 16-012 sections at a = 4 • angle of attack found in the literature.
It is important to note that the proposed optimization scheme is not subjected to any constraints, and therefore it is essential to verify during the optimization process that the attached cavity remains compatible with the physical assumption of the free-streamline problem, i.e., the attached cavity boundary does not intersect with the original hydrofoil boundary. The only intersections allowed are the detachment and reattachment points of the cavity. This issue can be resolved with the use of penalty functions to ensure that the attached cavity geometry is an admissible solution and does not intersect the original hydrofoil geometry between the detachment and reattachment points. Algorithmic refinements in future versions could include improving the convergence rate of the adjoint optimizer and enriching the cost functional with more terms in order to "flatten" the predicted pressure profile.
Overall, this method has been shown to predict well leading-edge cavitation; however, this methodology can also be extended to treat mid-chord cavitation. The present method is also suitable for the prediction of 3D sheet cavitation, strip-wise, when the cavity length chord-wise is a given quantity. Future work is planned toward a more systematic comparison of the present method with experimental data and other methods, as well as the investigation of the effects of camber on cavity shape and volume. Treatment of the 'direct problem', where the cavitation number is known, and the cavity shape and length are to be determined upon solution of the problem, is a challenging variation left for future work. In addition, the present model within the limitations of the ideal flow assumptions could be directly extended to treat the problem of cavitating hydrofoils operating beneath the free surface.