2.1 Blade Element Momentum Theory (BEMT) and tip loss model

The aerodynamic coefficients of the rotor in hover were predicted using BEMT. In this approach, momentum theory was employed to determine the induced inflow, while blade element theory was used to evaluate the two-dimensional aerodynamic characteristics at each radial blade section. The rotor disk was discretised into a finite number of radial elements, and the local aerodynamic quantities were computed iteratively for each element. The overall thrust and power coefficients were then obtained by integrating the sectional aerodynamic forces along the blade span [Reference Leishman2, Reference Johnson3, Reference Prouty11].

2.1.1 Induced inflow and Prandtl tip loss interaction

The induced inflow ratio ( $\lambda$ ), which varies along the rotor radius, is coupled with the Prandtl tip-loss factor ( $F$ ) which accounts for the reduction of vortex-induced effects near the blade tip [Reference Leishman2]. This condition ensures that the lift generated by the rotor blades approaches zero toward the tip region and is typically solved iteratively. In the present formulation, the momentum balance is applied locally at each annular radial station of the rotor disk. Accordingly, the induced inflow ratio is evaluated as a radially varying local quantity $\lambda (r)$ , rather than by prescribing a single uniform inflow value over the rotor disk. The non-dimensional induced inflow ratio at each radial station is given in Equation (1):

(1) \begin{equation}\lambda \!\left(r\right)=\frac{\sigma \!\left(r\right)C_{{l_{\alpha ,M}}}\left(r\right)}{16F\!\left(r\right)}\left(\sqrt{1+\frac{32F\left(r\right)}{\sigma \!\left(r\right)C_{{l_{\alpha ,M}}}\left(r\right)}\theta \left(r\right)r}-1\right)\end{equation}

The tip loss factor $F(r)$ in this equation is calculated using the Prandtl model, as given in Equations (2) and (3), depending on the radial position and the instantaneous induced inflow ratio:

(2) \begin{equation}f\!\left(r\right)=\frac{N_{b}}{2}\left(\frac{1-r}{\lambda \left(r\right)}\right)\end{equation}

(3) \begin{equation}F\!\left(r\right)=\frac{2}{\pi }\arccos \!\left(\exp \!\left(-f\left(r\right)\right)\right)\end{equation}

In this approach, Equations (1) and (3) form a coupled local system through the induced inflow ratio $\lambda (r)$ and the Prandtl tip-loss factor, $F(r)$ . These two quantities are updated successively at each radial station until the relative change between two successive values of $\lambda$ satisfies the convergence tolerance defined in Equation (4). The converged local solution is then used for the corresponding blade element before the calculation proceeds to the next radial station.

(4) \begin{equation}\frac{\left| \lambda _{i+1}-\lambda _{i}\right| }{\lambda _{i}}\lt 10^{-5}\end{equation}

The induced inflow ratio distribution obtained from calculations, varying along the blades as a function of the non-dimensional radius (r/R), is shown in Figure 2.

Figure 2.

Radial distribution of induced inflow (8 degrees).

2.1.2 Kinematic relations, compressibility effects, and viscous separation constraints

In BEMT analysis, the local angle-of-attack ( $\alpha$ ) on the blade element is defined as the difference between the local blade pitch angle ( $\theta$ ) and the inflow angle ( $\phi$ ) induced by rotor inflow [Reference Leishman2, Reference Johnson3]. Under the small-angle assumption, the inflow angle and the effective angle-of-attack are given by Equations (5) and (6), respectively:

(5) \begin{equation}\phi \left(r\right)=\frac{\lambda \left(r\right)}{r}\end{equation}

(6) \begin{equation}\alpha \left(r\right)=\theta \left(r\right)-\phi \left(r\right)\end{equation}

Due to the design characteristics of the XV-15 aircraft, the high rotor tip Mach number ( $M_{tip}=0.69$ ) makes compressibility significant at the outer radial stations. To model this variation in lift characteristics, the two-dimensional linear lift curve slope ( $C_{{l_{\alpha }}}$ ) was corrected using Equation (7) as a function of the local Mach number $M(r)$ in accordance with the Prandtl-Glauert approach [Reference Anderson12].

(7) \begin{equation}C_{{l_{\alpha ,M}}}\left(r\right)=\frac{C_{{l_{\alpha }}}}{\sqrt{1-M\left(r\right)^{2}}}\end{equation}

Using the corrected aerodynamic slope, the local lift coefficient ( $C_{l}$ ) at each radial station was calculated using Equation (8), with reference to the zero-lift angle-of-attack ( $\alpha _{ZLL}$ ):

(8) \begin{equation}C_{l}\left(r\right)=C_{{l_{\alpha ,M}}}\left(r\right)\left(\alpha \left(r\right)-\alpha _{ZLL}\right)\end{equation}

After the local lift coefficient is obtained, the profile drag coefficient of the blade element is evaluated. At high angles of attack, a rapid increase in drag is observed due to flow separation over the aerofoil [Reference Anderson12]. The high geometric twist angles near the root region of the XV-15 rotor increase the likelihood of viscous separation in hover conditions. To numerically represent this aerodynamic behaviour, it is assumed that the aerofoil enters the stall regime if the local lift coefficient exceeds the limits of the linear attached-flow region ( $-0.8\leq C_{l}\leq 1.2$ ) [Reference Leishman2]. Accordingly, the local drag coefficient was calculated using an empirical viscous separation model formulated by Equation (9):

(9) \begin{equation}C_{d}\left(r\right)=\begin{cases} C_{{d_{0}}}+k_{drag}C_{l}\left(r\right)^{2}, &-0.8\leq C_{l}\left(r\right)\leq 1.2\\ C_{{d_{0}}}+k_{drag}C_{l}\left(r\right)^{2}+\Delta C_{d,\textit{stall}} &\quad\qquad \textit{other}\,case \end{cases}\end{equation}

In this equation, $C_{d0}$ denotes the profile drag at zero lift, $k_{drag}$ is the drag rise coefficient, and ${\unicode{x0394}} C_{d,\textit{stall}}$ represents the abrupt increase in viscous drag due to flow separation.

2.1.3 Numerical trim and performance integration

In the final stage of the BEMT analysis, the differential forces calculated on each radial element ( $r_{k}$ ) were numerically integrated along the blade to determine the total rotor performance. The total thrust coefficient and power coefficient of the rotor were numerically summed over the step interval ${\unicode{x0394}} (r)$ between the radial stations using Equations (10) and (11), respectively [Reference Leishman2].

(10) \begin{equation}C_{T}={\sum }_{k=1}^{N}\frac{1}{2}\sigma \left(r_{k}\right)C_{l}\left(r_{k}\right){r}_{k}^{2}{\unicode{x0394}} r\end{equation}

(11) \begin{equation}C_{P}={\sum }_{k=1}^{N}\frac{1}{2}\sigma \left(r_{k}\right)\left(C_{l}\left(r_{k}\right)\lambda \left(r_{k}\right)+C_{d}\left(r_{k}\right)r_{k}\right){r}_{k}^{2}{\unicode{x0394}} r\end{equation}

In this study, to ensure sufficient numerical accuracy in the first stage, the radial domain was discretised into $N=40$ equal elements $({\unicode{x0394}} r=1/40)$ .

For the nominal DGW hover condition of the XV-15, the target thrust coefficient was set to ( $C_{T,\textit{target}}=0.00939$ ). This value corresponds to an isolated-rotor approximation in which the 13,000-lbf DGW [Reference Marr18] is equally shared between the two rotors, while wing/airframe download and rotor–airframe interactions are neglected. For each candidate twist distribution, the collective-pitch setting, denoted by $\theta _{root}$ in the trim formulation, was adjusted to satisfy this thrust requirement. For this purpose, the initial collective pitch angle was estimated using momentum theory (Equation (12)), and a trim loop was implemented and updated at each iteration. The residual error between the computed and target thrust $(\epsilon _{trim}=C_{T,\textit{target}}-C_{T})$ was evaluated, and the collective-pitch setting was updated using Equation (13):

(12) \begin{equation}{\unicode{x03B8}} _{{0_{init}}}=\frac{6C_{T}}{{\unicode{x03C3}} a}+\frac{3}{2}\sqrt{\frac{C_{T}}{2}}-\frac{3}{4}{\unicode{x03B8}} _{tw}\end{equation}

(13) \begin{equation}\theta _{root,i+1}=\theta _{root,i}+K_{trim}(C_{T,\textit{target}}-C_{T})\end{equation}

In Equation (13), $K_{trim}$ represents the convergence coefficient ( $K_{trim}=2.5$ ). The trim loop was continued until the error fell below the tolerance of $10^{-5}$ or the maximum limit of 200 iterations was reached. This approach ensures that the power consumptions of all investigated aerodynamic designs are evaluated in a consistent and comparable manner under conditions where they perform the same lifting task. To verify the physical consistency of the model and boundary conditions, the radial thrust distribution obtained for $\theta _{root}=8^{\circ }$ is presented in Figure 3.

Figure 3.

Radial distribution of thrust coefficient (8 degrees).

The distribution profile shows that the aerodynamic load is concentrated in the outer region due to the increasing dynamic pressure and decreases at the tip boundary $(r\rightarrow R)$ in accordance with physical expectations due to the Prandtl tip loss effect [Reference Leishman2].

2.2 Numerical solution and grid independence

In the developed BEMT model, the total thrust and power coefficients were obtained by numerically integrating the differential aerodynamic forces on the blade. For this purpose, the rotor disk was discretised into a finite number of radial elements. A low number of elements increases discretisation error, whereas excessively high numbers significantly increase computational cost, particularly in optimisation loops. Therefore, a grid independence study was performed to determine an appropriate number of radial elements.

The grid independence analysis performed for the hover condition $(\theta _{root}=5^{\circ })$ of the XV-15 rotor is presented in Figure 4. In this study, the number of radial elements was varied between $N=10$ and $N=100$ . The thrust values predicted by the model were compared with a high-resolution reference solution ( $N=150)$ , and the relative discretisation error was calculated.

Figure 4.

Grid independence and computational cost analysis.

The results show discretisation error is significant for $N\lt 30$ range. When the number of elements reaches N = 40, the relative error decreases to approximately 0.3 %. Beyond this value, further refinement yields only marginal improvements in the solution. Increasing the number of elements beyond N = 40 does not provide a significant gain in accuracy, while it increases computational cost. Therefore, N = 40 was selected for all subsequent numerical analyses and optimisation studies.

2.3 Model validation: comparison with NASA experimental data

The XV-15 tiltrotor experimental data published by NASA’s Ames Research Center were used as a reference for the validation of the numerical model. The comparison was performed using the Run 14 data set, which consists of isolated rotor hover tests in which rotor–fuselage interactions are neglected. Consistent with the experimental conditions, the analysis was conducted at standard sea-level conditions where the tip Mach number is $M_{tip}=0.69$ [Reference Harris10].

To evaluate the response of the numerical model under different aerodynamic loading conditions, the root collective angle was varied in the range of $\,1^{\circ }$ to $11^{\circ }$ . The resulting thrust coefficient, power coefficient and figure of merit (FOM) values were compared with the experimental data and are presented in Table 2.

Table 2.

BEMT model validation (comparison with NASA XV-15)

Table 2. Long description

A table comparing BEMT model validation with NASA XV-15 data. The table has 7 rows and 10 columns. Column headers are CT, CT, CP, CP error, FM, FM, Coll. CT (Ref) error (%), CP (Ref) (%), and FM (Ref) error (%). Row 1: 0.00482, 0.00473, 0.00039, 0.00038, 0.6076, 0.6025, 1.8, 1.9, 0.8. Row 2: 0.00624, 0.00641, 0.00047, 0.00052, 0.7444, 0.6957, 2.6, 10.7, 7.0. Row 3: 0.00750, 0.00815, 0.00062, 0.00070, 0.7461, 0.7398, 7.9, 12.4, 0.8. Row 4: 0.00987, 0.01096, 0.00093, 0.00108, 0.7494, 0.7540, 9.9, 14.0, 0.6. Row 5: 0.01168, 0.01301, 0.00119, 0.00143, 0.7501, 0.7357, 10.3, 16.6, 2.0. Row 6: 0.01260, 0.01398, 0.00134, 0.00161, 0.7465, 0.7236, 9.9, 17.1, 3.2.

The results in Table 2 show that the numerical model captures the overall aerodynamic trend. The error levels remain limited at low collective angles ( $\theta _{root}\leq 3^{\circ }$ ). In contrast, under medium and high lift conditions ( $5^{\circ }$ – $11^{\circ }$ ), the BEMT model underpredicts the thrust by approximately 7–10% compared to the experimental data. The power coefficient is also underpredicted by approximately 12–17%.

The nominal DGW hover target used in the optimisation ( $C_{T,\textit{target}}=0.00939$ ) lies within the loading range examined in the validation and is close to the BEMT result obtained at 8°collective ( $C_{T}=0.00987$ ).

These deviations are primarily attributed to the two-dimensional aerodynamic formulation of the BEMT approach. In particular, three-dimensional flow separation at the rotor root region, tip vortex effects and increased induced drag at high angles of attack cannot be fully captured. As a result, the model predicts lower $C_{P}$ and $C_{T}$ values at medium and high collective angles relative to the experiments.

When the FOM is considered separately, the discrepancy between the numerical model and the experimental data is smaller compared to the $C_{P}$ and $C_{T}$ results.

The thrust and power coefficient graphs shown in Figure 5 indicate a systematic underprediction tendency, particularly at medium and high collective angles. In contrast, the FM curve follows the experimental trend more closely over the same range.

Figure 5.

Thrust, power and FOM validation plots.

Figure 5. Long description

Panel A: A line graph shows the thrust coefficient (C_T) on the vertical axis, measured in units of 10^-3, against the collective pitch angle on the horizontal axis, measured in degrees. The graph includes two data series: one from the BEMT Model represented by a blue line and another from NASA XV-15 (Experimental) represented by orange squares. The thrust coefficient increases with the collective pitch angle for both data series. Panel B: A line graph displays the power coefficient (C_P) on the vertical axis, measured in units of 10^-3, against the collective pitch angle on the horizontal axis, measured in degrees. The graph includes two data series: one from the BEMT Model represented by a blue line and another from NASA XV-15 (Experimental) represented by orange squares. The power coefficient also increases with the collective pitch angle for both data series. Panel C: A line graph illustrates the figure of merit (FM) on the vertical axis against the collective pitch angle on the horizontal axis, measured in degrees. The graph includes two data series: one from the BEMT Model represented by a blue line and another from NASA XV-15 (Experimental) represented by orange squares. The figure of merit initially increases and then stabilizes as the collective pitch angle increases.

Due to the relationship $FM\propto \frac{{C}_{T}^{1.5}}{C_{P}}$ , the error trends in $C_{T}$ and $C_{P}$ are partially compensated in the FOM calculation. As a result, the deviation in FOM values remains more limited compared to those of the thrust and power coefficients. Within the medium-to-high loading range examined in the validation ${5^{\circ }-11^{\circ }}$ , the numerical FOM predictions deviate from the experimental data by approximately 0.6–3.2%.

The $C_{T}-C_{P}$ relationship shown in Figure 6 indicates that the model captures the general trend between thrust and power coefficients. Although a relative underprediction is observed in both $C_{T}$ and $C_{P}$ at medium and high collective angles, the effect of these deviations remains limited in the FOM calculation. These validation results support the use of the BEMT solver for comparative design evaluations within the present low-order framework, while the systematic underprediction observed at medium and high loading remains an acknowledged model-form limitation.

Figure 6.

Thrust and power coefficient comparison analysis.

3.1 Design variables: non-linear twist polynomial

In conventional helicopter blades, linear twist distributions are widely used because of their geometric simplicity; however, they may not provide an aerodynamically favourable twist distribution over the entire blade span under varying local dynamic-pressure conditions [Reference Leishman2]. In the XV-15 aircraft, however, to ensure that blade sections operate at optimal angles of attack corresponding to the varying relative velocities along the radial direction, the blade twist geometry is represented using a fifth-degree polynomial as a function of the non-dimensional radial position $r=r_{\text{local}}/R$ [Reference Harris10].

(14) \begin{equation}{\unicode{x03B8}} _{\text{twist}}\left(r\right)=p_{5}r^{5}+p_{4}r^{4}+p_{3}r^{3}+p_{2}r^{2}+p_{1}r+p_{0}\end{equation}

In Equation (14), the independent decision vector is defined as $\boldsymbol{X}=[\,p_{5},p_{4},p_{3},p_{2},p_{1}]$ , while $p_{0}$ is determined from the reference condition ${\unicode{x03B8}} _{\text{twist}}(0.75)=0$ . Since $p_{0}$ is the constant term of the polynomial, it does not affect the spanwise shape or gradient of the twist distribution, but only shifts the entire distribution uniformly. As the overall pitch level is independently adjusted through the collective-pitch trim procedure, treating $p_{0}$ as an additional independent design variable would introduce a redundant degree of freedom (DoF). Accordingly, $p_{0}$ is evaluated as

(15) \begin{equation}p_{0}=-\left[p_{5}\left(0.75\right)^{5}+p_{4}\left(0.75\right)^{4}+p_{3}\left(0.75\right)^{3}+p_{2}\left(0.75\right)^{2}+p_{1}\left(0.75\right)\right]\end{equation}

The lower and upper bounds of the five independent polynomial coefficients were defined by applying ±5% variations about their corresponding reference values, with the bound ordering adjusted according to the sign of each coefficient.

3.2 Objective function: power minimisation under constant thrust

The fundamental indicator of hover performance in helicopter and tiltrotor platforms is the shaft power required to meet a specified takeoff weight. Accordingly, the objective function of the proposed optimisation algorithm is defined as the minimisation of the rotor aerodynamic power coefficient $C_{P}$ .

Each candidate blade twist geometry generated from the decision variables is evaluated through the aerodynamic equilibrium condition within the BEMT solver prior to power computation. To ensure that the target thrust coefficient corresponding to the nominal DGW hover condition of the XV-15 rotor ( $C_{T,\text{target}}=0.00939$ ) is satisfied, the blade root pitch angle is iteratively updated until the thrust error is minimised. When the rotor disk satisfies the target thrust condition ( $\mid C_{T}-C_{T,\text{target}}\mid \leq 1\times 10^{-5}$ ), the corresponding power coefficient is computed, and the fitness function $J(X)$ is evaluated as follows:

(16) \begin{equation}J\left(X\right)=C_{P}\left(X\right)={\sum }_{i=1}^{N}\frac{1}{2}{\unicode{x03C3}} _{i}\left(C_{l,i}{\unicode{x03BB}} _{i}+C_{d,i}r_{i}\right){r}_{i}^{2}{\unicode{x0394}} r\end{equation}

In Equation (16), $N$ denotes the number of blade elements along the spanwise direction [Reference Johnson3]. The optimisation process seeks the parameter set that minimises the cost function defined in Equation (16). Minimising the power coefficient under a constant thrust condition is equivalent to maximising the aerodynamic efficiency of the rotor.

3.3 Artificial bee colony (ABC) algorithm setup

In the aerodynamic optimisation of the rotor blade, the ABC algorithm, a population-based metaheuristic method, was employed to handle the multidimensional search process in the design space [Reference Karaboga8, Reference Karaboga and Basturk9]. The algorithm maintains candidate designs within the prescribed aerodynamic and geometric feasibility criteria while improving the aerodynamic performance of the XV-15 rotor and is based on the collective behaviour of employed, onlooker and scout bees.

The search process is initialised with a population in which each employed bee represents a randomly generated twist polynomial. During the iterative loop, each bee generates a new candidate solution in the neighbourhood of its current position. This exploration in the continuous design space is governed by Equation (17):

(17) \begin{equation}v_{i,\,j}=x_{i,\,j}+{\unicode{x03D5}} _{i,\,j}\left(x_{i,\,j}-x_{k,\,j}\right)\end{equation}

Here, $k$ denotes a randomly selected candidate solution from the population ( $k\neq i$ ), while $\phi _{i,\,j}$ is a random coefficient in the range [−1, 1] that controls the magnitude of the search step. Each newly generated blade geometry is evaluated using the BEMT solver, and the resulting aerodynamic power coefficient is transformed into an inversely proportional fitness value for the selection mechanism of the algorithm:

(18) \begin{equation}{f\kern-1pt it}_{i}=1/\left(1+C_{P,i}\right)\end{equation}

After the employed bees share the fitness values within the colony, the onlooker bees select candidate solutions based on a probabilistic selection strategy $P_{i}$ .

(19) \begin{equation}P_{i}=\frac{{f\kern-1pt it}_{i}}{{\sum }_{n=1}^{SN}\,{f\kern-1pt it}_{n}}\end{equation}

According to Equation (19), the search process is biased toward regions with a higher probability of yielding lower power coefficients. The scout bee mechanism, on the other hand, supports the global exploration capability of the algorithm by abandoning solutions that cannot be improved within the prescribed scout limit of 300 trials.

Because of the stochastic nature of the ABC algorithm, the single-objective optimisation was performed through 30 independent runs. Each run was initialised using a different pseudo-random seed, with a colony size of 50 and a maximum budget of 10,000 number of function evaluations (NFE) per run. The solution yielding the lowest final power coefficient among the feasible solutions obtained from these independent runs was selected as the best-found solution for the subsequent aerodynamic evaluations. The stochastic repeatability of the optimisation and the consistency of the resulting twist distributions are further examined in Section 4.4.1.

The components of the search algorithm constitute the underlying framework of the multi-objective Pareto optimisation, which is presented in detail in Section 3.4.

3.4 Multi-objective optimisation (MOABC) and Pareto approach

In the second optimisation stage, total root-to-tip blade twist, ${\unicode{x0394}} \theta =| \theta _{\text{root}}-\theta _{\mathrm{tip}}|$ , was introduced as a second objective in addition to the aerodynamic power coefficient, $C_{P}$ . This geometric quantity provides a scalar representation of the overall blade-twist level, enabling a systematic examination of the trade-off between aerodynamic performance and alternative designs with lower total twist. The manufacturing-related physical basis for this choice is supported by composite-blade studies showing that root-to-tip twist and blade geometry are associated with fibre-path development during fabrication and with mould configuration [Reference Ho-Nguyen-Tan, Hong, Jaiswal, Oh, Bae, Yu, Kim and Yoon32, Reference Torpan and Zaharia34]. Nevertheless, manufacturing behaviour also depends on local variables such as surface curvature, fibre paths, steering radius and process conditions [Reference Bakhshi and Hojjati33]; therefore, ${\unicode{x0394}} \theta$ is used here as a geometric objective representing total blade twist rather than as a direct measure of manufacturability.

Within this framework, two design objectives are considered simultaneously. The first objective function, $f_{1}$ , is the minimisation of the aerodynamic power coefficient, $C_{P}$ , whereas the second objective function, $f_{2}$ , is the minimisation of the total root-to-tip blade twist variation, ${\unicode{x0394}} \theta$ . These conflicting objectives are transformed into a scalar cost function using the weighted sum method [Reference Gunantara13, Reference Marler and Arora14]:

(20) \begin{equation}J_{\textit{total}}\left(X\right)=w_{1}\cdot \overline{f}_{\!1}\left(X\right)+w_{2}\cdot \overline{f}_{\!2}\left(X\right)\end{equation}

Here, $w_{1}$ and $w_{2}$ denote the weight coefficients that vary according to design priorities $(w_{1}+w_{2}=1)$ , and $\overline{f}_{\!n}$ represents the normalisation values that enable objectives with different orders of magnitude to be summed together. The normalisation process was performed by squeezing both objectives between 0 and 1:

(21) \begin{equation}\overline{f}_{\!n}=\frac{f_{n}-f_{n,min}}{f_{n,max}-f_{n,min}}\end{equation}

The algorithm was run by systematically varying the weighting coefficients, resulting in a Pareto front rather than a single solution. This front represents the trade-off between the aerodynamic power coefficient and total blade twist and defines the non-dominated solution set with respect to $C_{P}$ and ${\unicode{x0394}} \theta$ [Reference Gunantara13, Reference Marler and Arora14].

Through this approach, the design alternatives presented in Section 4 are based on a consistent analytical framework that supports engineering decision-making.