An Efficient Parametric Modeling, Evaluation and Optimization Strategy for Aerodynamic Configuration Design of eVTOL Aircraft

Abstract

Traditional aerodynamic design faces significant limitations in modeling and computational efficiency during conceptual design stage. A phased collaborative aerodynamic design strategy for eVTOL aircraft was established, by combining the OpenVSP platform for rapid parametric modeling and evaluation, a Kriging surrogate framework with an improved differential evolution algorithm for optimization, and the SUAVE platform for propeller reverse design. In the wing-body (WB) optimization phase, 23 configuration parameters such as the wing shape and location were adjusted. The aerodynamic evaluation was conducted using the Vortex Lattice Method (VLM) in OpenVSP, resulting in a 9.3% increase in the lift-to-drag ratio (L/D). During the wing-body-propeller (WBP) coupling optimization phase, the Actuator Disk Theory (ADT) was incorporated into WB model to quantify the slipstream effects. After optimizing the key geometric parameters such as disk diameter and location, the comprehensive propulsion efficiency and lift-to-drag ratio (η·/L/D) was increased by 14%. Relative performance parameters were then transferred to SUAVE to reconstruct the propeller based on the Betz-BEM theory. The RANS high-fidelity verification of the optimized WBP model shows high consistency in the trends of lift coefficient C_l and L/D calculated by VLM, with the propeller thrust error 5.2%, and the C_l error 9.7%, which confirms the engineering reliability and efficiency of the proposed strategy.

Introduction

With urbanization accelerating, traditional ground transportation systems face severe traffic congestion, with peak-hour vehicle speeds below 15 km/h in megacities [1]. The Urban Air Mobility (UAM) systems have emerged as a breakthrough solution, with Electric Vertical Take-Off and Landing (eVTOL) aircraft becoming a pivotal technology due to its vertical takeoff/landing capability, high cruising efficiency, and environmental benefits [2–5].

Aircraft modeling is the primary link in aircraft design. Although traditional CAD tools (such as CATIA) are widely used, their limitations become prominent in the conceptual design stage. Ronzheimer et al. [6] developed a design table-driven parametric framework using CATIA V5, enabling high-precision modeling and CFD-compatible optimization of propeller blades and blended-wing-body (BWB) configurations. Bhise et al. [7] proposed a parametric CATIA modeling approach for wings of the short range and narrow body aircraft, integrating Ansys Fluent aerodynamic analysis with Python scripting to generate NURBS curves. However, these methodologies rely on manual hierarchical feature tree construction, resulting in disconnected parameterization logic from engineering semantics. Geometric inconsistencies frequently arise during parameter adjustments, and iterative efficiency remains inadequate for large-scale design exploration. Regarding deformation-based methods, Lyu et al. [8] employed Free-Form Deformation (FFD) for aerodynamic optimization of BWB configurations, quantifying design variable impacts on performance metrics. Li et al. [9] combined FFD with a discrete adjoint method using Reynolds-Averaged Navier-Stokes (RANS) simulations, achieving a 5.7% reduction in total drag for supersonic airliner configurations. Nevertheless, FFD methods suffer from the “curse of dimensionality” due to excessive control points, limiting their applicability to localized geometry adjustments. Their dependence on initial baseline models prevents autonomous topology generation, thereby constraining innovative design possibilities.

As the core power component of the eVTOL aircraft, the aerodynamic characteristic analysis of the propeller is the basis for design optimization. For performance evaluation of an isolated propeller, high-fidelity Computational Fluid Dynamics (CFD) methods based on the RANS equations are widely utilized, such as the Multiple Reference Frame (MRF)-RANS and Unsteady RANS (URANS). These methods are capable of accurately capturing three-dimensional flow and viscous effects. For instance, Bryant et al. [10] utilized the RANS method to calculate the thrust of a small-scale propeller, but it took 14–15 h (160 threads) under a single working condition. Liu et al. [11] compared the accuracy differences between MRF-RANS and the Blade Element Momentum Theory (BEMT), and found that while RANS exhibited lower errors (with propulsive efficiency errors below 10%), it required 900 CPU hours under a single working condition. However, in the actual flight conditions of eVTOL aircraft, the strong coupling effects between the propeller and the wing significantly influence the overall aerodynamic performance. Consequently, analyzing only the isolated propeller is insufficient to meet design requirements.

For the wing-propeller aerodynamic interaction analysis, it is also necessary to balance the efficiency and accuracy of high-fidelity methods. Sinnige et al. [12] employed the Unsteady Reynolds-Averaged Navier-Stokes (URANS) approach to analyze the dynamic stall characteristics of the NASA LA-8 tilt-configuration aircraft during cruise and transition flight states. The error compared with experimental data was less than 8%. Li Peng et al. [13] explored the influence of tilting rotors on the wing and fuselage drag in hovering and cruising states through the RANS-based overset grid technique. Xu et al. [14] analyzed the influence of the lower slipstream of the upper propeller of the coaxial reverse-rotating propeller configuration eVTOL aircraft on the pull performance of the lower propeller based on MRF-RANS, and quantified the pull performance loss (about 33%). Although the RANS method can analyze the nonlinear influence of slipstream on the pressure distribution over the wing [15–18], its single-condition calculation cost is as high as several hundred CPU hours, which seriously restricts the optimization efficiency.

In the field of aerodynamic shape optimization, although the research based on surrogate models or discrete adjoint gradient methods has achieved the coupling optimization of propellers and wings [19–21], aerodynamic analyses all rely on high-fidelity calculation methods such as RASN. For example, Yildirim [22] utilized a high-fidelity computational method based on the OpenMDAO/MPhys framework for the optimization of an eVTOL wing-propeller aerodynamic interaction, which achieved an 18.3% reduction in power under all constraint conditions. Chauhan [23] conducted aerodynamic simulations of a distributed electric propulsion eVTOL using the ADT-RANS method, which combines the Actuator Disk Theory (ADT) with Reynolds-Averaged Navier-Stokes (RANS) equations. A gradient-based optimization algorithm was then employed to minimize cruise drag. However, the design variables were limited to the propeller’s rotation direction, the wing’s twist, and airfoil shape. Pedreiro [24] referred to the Embraer EMB-120 and optimized the twist distribution and airfoil shape of the wing behind the propeller using RANS-ADT with six design variables. Eventually, the drag was reduced by three counts under cruising conditions. Although the high-fidelity method based on the RANS equations can accurately simulate the viscous interference effects between the propeller slipstream and the wing, the significantly high computational cost hinders multivariable optimization of eVTOL. To reduce complexity, while there are existing studies have combined low-fidelity vortex lattice methods or BEMT with surrogate models for aerodynamic shape optimization of eVTOL, the optimizations of the wing and propeller are performed separately [25, 26].

To balance the contradiction between efficiency and accuracy in aircraft design and aerodynamic optimization, the open-source tools based on NASA’s OpenVSP offer a new approach. Its parametric component library and script-driven modeling framework enable the construction of eVTOL models within minutes through adjusting predefined geometric variables, eliminating the need for manual feature reconstruction at each design iteration. Its Unsteady Vortex Lattice Method (UVLM) enables rapid performance evaluation of isolated propellers. For instance, Sun Zongyan et al. [27] adopted UVLM to analyze the distributed propulsion configuration eVTOL, achieving lift coefficient C_l and drag coefficient Cd errors less than 5% and 4.6% within the linear region. Yu et al. [28] utilized OpenVSP’s UVLM and the steady RANS method based on the k-ω SST turbulence model in Simcenter STAR-CCM+ software to conduct performance analyses on a single three-bladed propeller. The differences in thrust and torque calculations between the two methods were less than 10% and 5%. Sheridan et al. [29] calculated propeller thrust using both OpenVSP’ss UVLM and OVERFLOW’s URANS. The differences in thrust and torque calculations between the two methods were less than 10% and 2%. The UVLM method took only a few minutes under a single working condition. Wu et al. [30] developed the ADT-VLM (Actuator Disk Theory-Vortex Lattice Method) framework to enable efficient simulation of propeller-wing aerodynamic interactions. This approach achieves C_l and Cd deviations below 10% compared to RANS results. Yu et al. [28] validated ADT-VLM’s capability to capture three-bladed propeller slipstream characteristics, demonstrating C_l deviations below 6% relative to steady RANS simulations, while computational efficiency achieved less than 1 min under a single working condition.

However, existing studies have not yet fully utilized OpenVSP’s parametric modeling capabilities. There remains an urgent need to establish an integrated optimization framework combining ADT-VLM with surrogate modeling. Unlike conventional decoupled optimization strategies, the ADT-VLM methodology eliminates velocity field discrepancies between slipstream conditions and design points through coupled analysis. When integrated with Kriging surrogate models, this approach significantly reduces computational costs to achieve global propeller-wing coupling optimization.

Numerical Simulation Methods and Aerodynamic Validation

Vortex Lattice Method Theory

The Vortex Lattice Method (VLM) is based on the potential flow theory, which discretizes the aircraft surface into vortex lattices to simulate the surrounding flow field. It is suitable for rapid aerodynamic performance evaluation under small angle of attack and low speed conditions, offering significant efficiency advantages during conceptual design phases. OpenVSP implements VLM calculations through VSPAERO tool, which can output key parameters such as C_l and C_d.

Based on the lifting surface theory [31], the classical VLM simplifies the wing into discrete thin lifting surfaces with horseshoe vortices deployed. Each panel contains a horseshoe vortex composed of bound vortex lines and free vortex lines, with circulation strength Γ. A control point P is established at the midpoint of the 3/4 chord line of each panel (as shown in Figure 1a). The induced velocity from vortex segment ab at control point P is calculated using the Biot-Savart law [32], expressed as Equation 1:

FIGURE 1

In contrast, OpenVSP employs an improved VLM that replaces horseshoe vortices with a vortex ring model. Vortex rings are placed at panel centers, with their geometric center point P acting as the singularity (similar to the control point in horseshoe vortices). Wake vortices are represented by trailing vortex lines from trailing-edge wake points (Figure 1b). The Prandtl-Glauert compressibility correction factor is introduced to account for compressibility effects [33]. Induced velocities are computed using the Biot-Savart law, expressed as Equation 2:where represents the strength of the vortex ring, and K is the factor considering the influence of compressibility on the calculation: if , then K = 2; otherwise, K = 1, is defined by Equation 3

By incorporating these velocities into the Neumann boundary condition [34], the vortex ring strengths are solved iteratively. Aerodynamic loads are subsequently calculated using the Kutta-Joukowski theorem [35]. Some researchers contend that VLM implementation in VSPAERO constitutes an advanced lifting-surface method rather than a technically pure VLM [36].

Actuator Disk Theory

The Rankine-Froude momentum theory was extended to the actuator disk theory to estimate propeller thrust through simplified flow modeling [37]. VSPAERO implements an improved formulation combining Conway’s elliptical actuator disk model [38] and Johnson’s tangential velocity model [39], which enables the calculation of axial velocity, radial velocity, and the pressure jump across the disk surface.

According to [38], the axial velocity induced by a propeller of radius R at a point (r,z), where r is the radial coordinate ranging from 0 to R and z is the axial coordinate with the origin (0,0) (corresponding to the disk center) is described by Equations 4, 5:

The radial velocity is given by Equation 6, is defined by Equation 7, and the pressure jump across the disk surface is given by Equation 8.

In Equations 4–8, the remaining terms are defined as follows: the axial velocity at the center plane of the disk is given by Equation 9, and the induced velocity is given by Equation 10,where T is the propeller thrust (obtained by inputting the thrust coefficient C_T), and A is the projected area of the propeller disk. Finally, the parameter α (not to be confused with the angle of attack) is defined by Equation 11.

The tangential velocity is modeled considering the power losses due to the viscous drag of the blade. Consequently, Johnson’s equation relies on the blade’s lift-to-drag ratio (C_l/C_d). However, VSPAERO, based on the assumptions outlined in literature [39], circumvents the need to provide blade aerodynamic data. The tangential velocity is given by Equation 12, where is defined by Equation 13.

The C_l/C_d is modeled by combining the blade hover thrust coefficient () and the hover power coefficient (), leading to Equations 14, 15.

Equation 14 defines the average C_l assuming the entire blade is in operation, with σ representing the blade solidity. The first term in Equation 15 is the profile drag power coefficient for a blade with a constant chord length and a constant drag coefficient, and the second term represents the hover induced power loss, which is corrected by an empirical factor κ [39]. According to [39], when considering a linear inflow distribution, κ takes a value of approximately 1.17, which is also the value adopted in the VSPAERO code.

If input parameters C_T, C_P and RPM are provided, and can be calculated. Subsequently, the C_d/C_l term in Equation 12 can be estimated through simple Equation recombination, and then the tangential velocity can be calculated.

Grid Convergence Analysis

The computing resource in this work is the Intel Xeon E5-2698v3 processor, with a main frequency of 2.30 GHz, 16 cores and 32 threads. The method verification and optimization of OpenVSP based on the VLM, and the high-precision verification of RANS both use 20 threads.

The convergence of the grid is analyzed before the verification of VLM method. Based on OpenVSP’s VLM grid control logic, it employs an X-Z planar symmetric modeling method for wings, and the wing surface grid is controlled by chordwise nodes Num_W and spanwise nodes Num_U. The grid quantities in the chord direction and span direction are (Num_W-1)/2 and Num_U-1 respectively. To ensure mathematical closure of grid topology, OpenVSP requires Num_W = 4n+1 for wings and Num_W = 8n+1 for propellers (n is a positive integer). Additionally, grid refinement can be applied at specific locations such as the wingtip, wingroot, leading edge, or trailing edge. When the value is between 0 and 1, the grid will be refined, and the smaller the value, the higher the grid refinement level. When the grid refinement value is set to 1, it indicates no grid refinement will be performed.

This section selects the classic symmetric airfoil NACA 0015 as a case study [40]. The computational model adopts a rectangular wing configuration with a chord length c = 1.7 feet (0.518 m), a semi-span length b/2 = 5.59 feet (1.704 m), and an aspect ratio AR = 6.45. The wingtip features a rounded transition design to simulate the geometric characteristics of an actual wind tunnel model (Figure 2), effectively eliminating tip effects caused by sharp edges. The numerical calculations strictly adhere to the experimental conditions: a freestream Mach number Ma = 0.13 ( m/s), Reynolds number Re = 1.5 × 10⁶, and air density ρ = 1.225 kg/m³.

FIGURE 2

The study of grid convergence mainly investigates the influences of three grid-related variables: total grid count, spanwise-to-chordwise grid ratio, and local grid refinement. For the spanwise to chordwise grid ratio, two configurations are analyzed: 2:1 and 1:1. Both ratios cover the refined grid division (the wing leading edge/tip/root refinement factor is 0.25, the trailing edge refinement factor is 1) and unrefined grid division (refinement factors uniformly are 1). A total of 26 grid configurations are generated, with grid counts ranging from 96 to 13,068. Typical grid distributions for refined and unrefined cases are illustrated in Figure 3. By comparing the convergence curves of the C_l and C_d for different grid parameters (as depicted in Figure 4), it was found that the 1:1 aspect ratio grid without refinement exhibited the optimal convergence characteristics. When the angle of attack is 10, the difference in C_l between the 6,144 grids and 13,616 grids cases was only 0.54%. For the 2:1 aspect ratio grid with refinement, a C_l fluctuation of 1.8% was observed at high grid refinements, indicating that overly refined grids near the leading edge might trigger numerical oscillations. Based on these findings, the unrefined grid with an aspect ratio of 1:1 and 97 spanwise cells × 65 chordwise cells (totaling 6,144 grid cells) is selected as the baseline grid configuration.

FIGURE 3

FIGURE 4

VLM Validation

Based on the above operating conditions and baseline grid configuration, the reliability of OpenVSP’s VLM is verified by calculating the C_l and C_d at different angles of attack in NACA 0015. Figure 5 shows the comparison of the calculation results of the VLM with experimental data in the 2°–14° angle of attack range. Within the 2°–12° angle of attack range, the computed C_l values exhibit a high degree of linear correlation with the experimental values (R² = 0.897), with a root mean square error (RMSE) of 0.125. At low angles of attack (α = 2°), the relative error is 14.0% (computed value 0.154 and experimental value 0.179). As the angle of attack increases, the error gradually rises. In the stall transition region (α = 14°), the predicted C_l is 1.055, which is 15.3% lower than the experimental value of 1.245. However, the trend maintains a monotonic increase, consistent with the experimental data and does not exhibit a stall inflection point. The phase error in predicting the stall angle of attack is less than 1° compared to the experimental data.

FIGURE 5

ADT-VLM Validation

The PROWIM model [41], which has comprehensive wind tunnel test data, was selected for validation of the ADT-VLM. The configuration of the model is illustrated in Figure 6. The main body of the model features a straight wing with an aspect ratio of 5.33, the airfoil is NACA 642-A015 and a semi-span length is 0.64 m (b/2). The propulsion system consists of a four-bladed metal propeller with diameter D = 0.236 m, and its power unit is positioned 0.3 m along the wing’s mean aerodynamic chord. The nacelle’s rotational axis is coplanar with the turntable’s reference plane, with the maximum output power P = 5.5 kW. The flow conditions as follows: Re = 8 × 10⁵, Ma = 0.15, ρ = 1.225 kg/m³, and propeller advance ratio J = 0.85. The unrefined grid with an aspect ratio of 1:1 and 66 spanwise cells × 49 chordwise cells (totaling 3,120 grid cells) is selected as the baseline grid configuration. The propeller slipstream interference effect is evaluated using ADT-VLM. A comparison between the computational results and the experimental data is presented in Figure 7.

FIGURE 6

FIGURE 7

Figure 7a shows that when α = 6°, the maximum C_l relative error between ADT-VLM and the test value is only 2.7%, which verifies the ability of this method to capture the slipstream enhancement effect. Notably, when α = 8°, the maximum relative error is merely 1.5%, indicating that it has high computational confidence at medium angle of attack. Figure 7b illustrates that when the propeller is on, the computed “C_d” appears negative because it includes the propeller’s thrust. By accounting for thrust correction and comparing the results with the experimental values, the validity of the method can still be verified. The results show that the maximum C_d error between ADT-VLM and the experimental value is less than 12.2%.

Comprehensive error analysis shows that VLM and ADT-VLM have good applicability in the concept design stage, and its single calculation time does not exceed 1 min, and its computing efficiency advantages are significant, laying the foundation for the subsequent optimization framework.

Rapid Parametric Modeling and Optimization Framework

OpenVSP-Based Parametric Modeling

In the process of parametric modeling and aerodynamic evaluation, the script-driven approach based on the OpenVSP platform achieves fully automated operations through predefined “.vspscript” files, as depicted in Figure 8. This script file embeds the geometric parameterization code and aerodynamic analysis code of each component of the aircraft. After Windows command line interface CMD calls the OpenVSP program, the script is executed to sequentially complete geometric modeling and aerodynamic analysis.

FIGURE 8

The fuselage modeling is rapidly accomplished using the Fuselage component to create a streamlined shape, with side and front views shown in Figure 9a. The fuselage is divided into four segments, and each segment has an axial position controlled by XLocPercent parameters. The cross-sectional shape adopts an elliptical profile (the aspect ratio defined by Ellipse_Width and Ellipse_Height), and the end closure method is optimized by CapUMaxOption parameter.

FIGURE 9

Wing modeling is implemented by Wing components, and its main parameters include Area, Span, Chord, Aspect Ratio, Taper Ratio, Sweep, Dihedral, Incidence and Twist. Figure 9b shows main wing is a straight wing with two segments of wings, which requires at least three wing shapes, three chord lengths and two spread lengths. The tail wing is a single-stage wing with a sweep angle and a dihedral angle. Its airfoil is defined by parameterization of the airfoil camber (Camber), The ratio of thickness to chord length (T/C) and maximum camber location (Camber Loc).

The propulsion system modeling includes two modes: the propeller mode (Blade Mode) and the actuator disk mode (Disk Mode), as respectively depicted in Figures 9c,d. In the propeller mode (Blade Mode), the Propeller component is utilized to define parameters such as the blade’s radial chord length, twist angle distribution, and airfoil distribution. This mode supports the import and export of Blade Element Momentum (BEM) files, facilitating collaborative design across multiple tools. This functionality was employed in the subsequent propeller reverse design.

The actuator disk mode (Disk Mode) simulates the aerodynamic interference effect of the propeller on the wing and fuselage and calculates the slipstream velocity increment by inputting parameters such as D, RPM, C_T, and C_P. This approach avoids the need for complex blade geometry modeling, thereby significantly reducing computational resource consumption.

For the aerodynamic analysis, the script file “.vspscript” defines parameters such as the angle of attack (Alpha), Mach number (Mach), and Reynolds number (ReCref) through the VSPAERO module, and specifies the VLM as the solution approach. After the calculation is completed, the geometric model file “.vsp” and the standardized aerodynamic report file “.polar” are automatically output. The file “.polar” contains key indicators such as C_l, C_d and induced drag coefficient C_di. Additionally, the model can be visualized by calling the “vspviewer” program via the CMD interface.

This process achieves full closed-loop automation of “geometric modeling- aerodynamic calculation-result output” through a single script file “.vspscript”. Each iteration takes less than 1 minute, which represents two orders of magnitude more efficiently than the traditional CFD method, providing reliable technical support for the rapid optimization of eVTOL aerodynamic configurations. The model file “.vsp” also supports direct export to format files such as STL, STP, and IGS. This method is used in the subsequent CFD high-fidelity verification.

Surrogate Optimization

This study constructed an optimization framework centered around the Kriging surrogate model in Figure 10a. This framework realizes efficient optimization of eVTOL aerodynamic configurations through a multi-layer collaboration mechanism.

FIGURE 10

In the initial sampling phase, an improved Latin hypercube sampling (LHS) method is employed to generate design space samples [42]. By integrating OpenVSP parametric modeling with VLM, aerodynamic evaluations and training dataset construction are completed within minutes. Subsequently, Gaussian process regression is utilized to establish a Kriging surrogate model [43], which characterizes the correlation between design variables and aerodynamic responses through nonlinear mapping relationships. This provides an efficient surrogate model for subsequent optimization processes.

The JADE-SPS-EIX algorithm [44] is employed as the optimization search algorithm, which introduces three innovations based on JADE [45]. Firstly, it adopts a dynamic parameter adaptation mechanism that adjusts the scaling factor (sf) and crossover probability (Cr) in real time using historical iteration data. Here, sf follows a Cauchy distribution, while Cr is updated based on a normal distribution derived from successful individuals, enabling the algorithm parameters to adaptively match the complex design space. Secondly, it constructs a Successful Parent Selection (SPS) strategy [46]. When the j individual in the population fails to produce a superior offspring in succession Q times, the algorithm forcibly selects a parent individual from the set of successful individuals (S) for mutation and crossover operations, thereby preventing the search from stagnation. Thirdly, it proposes an Eigenvector-based Crossover (EIX) operator [47, 48]. This operator constructs an eigenvector coordinate system through covariance matrix analysis and performs crossover operations within this system instead of the natural coordinate system. This approach allows the search direction to gradually approximate the feature space of the objective function, significantly enhancing the efficiency of solving coupled problems.

During the dynamic adding point and iterative updating phase, the Generalized Hybrid One-Two Stage Method (GHOTM) proposed in reference [44] is adopted to balance exploration and exploitation. This method employs a two-stage dynamic switching strategy: In the global exploration phase, based on the Generalized One Stage Method (GOM) [44], it comprehensively considers the predicted values of the surrogate model, prediction variances, and gradients of the objective function. It prioritizes regions with high uncertainty, low objective values, and smooth gradients as new sample points. When a potentially optimal region is detected, the algorithm automatically switches to the Minimize Prediction method (MP) [49] to perform local search. The method gradually narrows the search range centered around the current optimal solution for in-depth development. If the local search stagnates, the algorithm will automatically restart GOM for global exploration. To mitigate the impact of optimization boundary settings on the optimization results, a 10% adaptive boundary expansion method for the design space is employed.

This mechanism achieves an adaptive balance between global exploration breadth and local search depth by monitoring the optimization state in real time. The newly added sample points are validated twice by OpenVSP and VLM to update the surrogate model, forming an autonomous closed-loop iteration of “parametric modeling - surrogate optimization - convergence determination”. Through modular process design, the optimization framework cohesively integrates sampling, modeling, aerodynamic optimization and verification, ultimately achieving efficient optimization of complex aerodynamic configurations.

A Phased Optimization Strategy

This study adopts a phased collaborative aerodynamic design strategy in Figure 10b (“WB Optimization - WBP Coupling Optimization - Propeller Reverse Design - High-Fidelity CFD Validation”). This strategy systematically decouples complex aerodynamic coupling issues through hierarchical decomposition, significantly reducing the computational complexity of multivariable optimization. The first stage focuses on wing aerodynamic shape optimization, enhancing cruise efficiency through lift-to-drag ratio (L/D) maximization. The second stage introduces an actuator disk model to quantify propeller slipstream interference effects on the wing, employing a comprehensive performance target (η·L/D) as the optimization objective. The third stage performs reverse design to convert actuator disk parameters into propeller geometry model, with final validation conducted via high-fidelity CFD simulations to confirm the physical feasibility and engineering applicability of optimization results.

Optimization Results and Discussion

Wing-Body Optimization

To achieve the above goals, this study utilizes the Supernal S-A1 [50] (Tilting rotors fixed Lift propellers) model as a reference prototype and builds the baseline model as the optimization target in Figure 11a. Since this study focuses only on the cruise condition, the tilting propellers are modeled only in the propulsion state, ignoring the lift propeller interference in the non-working state. During the phase of WB optimization (Wing-Body), a simplified baseline model is adopted, which includes only the fuselage, main wing, and tail wing in Figure 11b. The geometric parameters of main wing are presented in Figure 9b.

FIGURE 11

The simplified baseline model references the actual flight conditions of the S-A1: flight altitude approximately H = 0.5 km, ρ = 1.16727 kg/m³, Ma = 0.169 (V = 57.18 m/s), Re = 4.39 × 10⁶, and α = 2°. The wing parametrization study is conducted with the optimization objective of maximum L/D. During the optimization process, the fuselage and tail configuration parameters (position and shape) remain fixed. Minor longitudinal adjustments (X-axis) of the wing position are permitted to maintain structural pitch moment characteristics, while minor vertical adjustments (Z-axis) are allowed to optimize ground effect performance. Aerodynamic configuration parameters of the wing include 23 design variables in Table 1. The pressure distribution is reconstructed and balances profile drag and induced drag through collaborative parameter adjustments. The mathematical model of aerodynamic optimization is established as Equation 16.

Among them, Cl_wing and Cl_wing,0 represent the lift coefficients before and after optimization, while S_wing and S_wing,0 denote the wing areas before and after optimization. Constraints include a lower limit for the C_l ≥ 0.638 and a fluctuation range for wing projected area (S_w ± 10%). Constraint violations are managed through a penalty function approach, with penalty factors set to 30 for area constraints and 600 for lift constraints to align with the objective function magnitude. LHS is employed to generate 30 initial samples, with a maximum evolutionary iteration (Maxit = 100) and a population size (NP = 100). Figure 14a shows efficient coordinated adjustment of aerodynamic configuration parameters is achieved through 11 rounds of adding point, including the coordinated iteration of global GOM and local MP. The sampling process takes 9 min, and the optimization process takes 3 h and 20 min, which indicates that the calculation time for a single working condition is less than 1 minute

In Table 1, optimization results indicate that the value of X_location 0.276 m exceeds the boundary value 0.25 by 10%, which is in line with the design space adaptive expansion mechanism. And the wingspan increased from 14.75m to 15.69 m (+6.4%), the mean chord length extended to 1.23 m (+2.1%), and a 0.603° sweep angle along with a −0.575° dihedral angle were introduced. For airfoil parameters, thickness ratio (T/C) of the outer-section airfoil (Airfoil 2) reduced from 0.11 to 0.10, while the camber value (Camber) increased from 0.04 to 0.06. These modifications reconstructed the leading-edge suction peak and trailing-edge pressure recovery gradient as shown in Figure 12. In Table 2, the optimized L/D improved from 28.9 to 31.6 (+9.3%), the C_l increased by 35.1% (from 0.638 to 0.862), and the wing area (S_w) only expanded by 7.5%, satisfying all constraints.

FIGURE 12

Wing-Body-Propeller Coupling Optimization

Based on the optimized WB configuration, ADT is employed to simplify the propeller into a disk. This approach is further used to conduct the propeller-wing coupling optimization through a new baseline WBP model (Wing-Body-Propeller), as shown in Figure 13. And the comprehensive propulsion efficiency and lift-to-drag ratio (η·/L/D) are taken as the optimization objectives, which aims to achieve the collaborative optimization of energy conversion efficiency and cruise aerodynamic performance.

FIGURE 13

With the optimized drag coefficient Cd = 0.0273 and drag Equation 17, the total cruise drag F is calculated as 935 N.

When total propulsion efficiency η_total = 0.75 is set, the thrust Equation 18 yields a total required total thrust T_total = 1246 N. It is subsequently distributed equally across four disks to obtain the single disk thrust T = 311.5 N.

During the optimization process, the shape and location of the wing remain constant. The design variables cover a total of 7 parameters, including the rotational speed (RPM), diameter (D), location and rotation angle of disks on the main wing and the tail wing, as shown in Table 3. With single disk thrust T = 311.5 N established as a fixed boundary condition, the C_T is derived using disk D and RPM as design variables. A dual-constraint mechanism governs C_P selection: first, an empirical feasible range C_P∈[0.04,0.08] is defined based on electric propulsion system efficiency characteristics, and then the theoretical boundary C_PEquation 19 is inversely deduced by the efficiency Equation.

The intersection of these constraints defines the physically realizable parameter space. During optimization iterations, the value of C_P needs to be secondary validation. When η > 0.9 or C_P is out of bounds, a forced penalty function of η = 0.5 is applied to ensure stable optimization. Consequently, the mathematical model of integrated propeller-wing aerodynamic optimization is established as Equation 20

Figure 14b shows the WBP optimization is achieved through 8 rounds of adding points, including the coordinated iteration of global GOM and local MP. The sampling process (30 samples) takes 7 min, and the optimization process (200 samples) takes 1 h and 50 min, which reveals that the calculation time for a single working condition does not exceed 1 min.

FIGURE 14

Table 3 shows that after optimization, the vertical position of the disk (Z_location) is moved up to 0.015 m (+0.14 m), the D is increased to 1.416 m (+4.1%), the rotation angle (Y_rotation) is adjusted to 3.736°, and the rotational speed is reduced to 2398 RPM. The optimized slipstream coverage pattern in Figure 15 reveals a lateral shift of the slipstream core toward the wing’s mean camber line. Concurrently, the leading-edge suction peak ΔCp_min undergoes pressure recovery from −7.474 to −3.692, the pressure gradient becomes more gradual, and the flow separation is effectively suppressed, which directly correspond to the 24.5% reduction of C_di in Table 4. This phenomenon conforms to the slipstream velocity correction Equation 21 in momentum theory.

FIGURE 15

The increase in diameter reduces both C_T and slipstream velocity . Despite a 4.1% reduction in C_l attributed to slipstream coverage adjustment, propulsive efficiency η demonstrates a 6.25% improvement. The η·L/D increases from 27.21 to 31.02 (+14%). The results show that by coordinating the location and diameter of the disk, the slipstream gain and aerodynamic interference loss can be effectively balanced, significantly improving the energy utilization efficiency of the eVTOL cruise phase.

Propeller Reverse Design

Based on the propeller-wing coupling optimization results, the SUAVE platform is used to realize the reverse reconstruction of disk parameters to three-dimensional propeller geometry based on Betz-BEM theory. The input parameters are shown in Table 5, and the radial distribution of the optimal propeller geometric characteristic parameters in this design state is shown in Figure 16a. The design propeller is evaluated by OpenVSP’s MRF-VLM method, and its pressure difference coefficient contour and vorticity contour are shown in Figures 16b,c.

FIGURE 16

The single propeller thrust calculated by MRF-VLM method is 301.0 N, which deviates 3.4% from SUAVE’s BEM -based design target value of 311.5 N. BEM provides rapid preliminary design and MRF-VLM enables three-dimensional flow verification, forming a design-validation closed loop. This validates the efficacy of the BEM-MRF hybrid design in eVTOL propeller inverse design.

High-Fidelity Validation

To systematically evaluate the physical rationality of the optimization design, STAR-CCM+ is used to carry out MRF-RANS numerical verification. For isolated propellers, a rotating and static domain partition grid is constructed in Figure 17a. Figure 17b shows that a total grid volume of 4.76 million is generated using non-structural grids, and grid refine is performed on the propeller blade tips to capture the boundary layer effect. Figure 17c shows its pressure coefficient cloud diagram.

FIGURE 17

The calculation results in Table 6 show that the thrust and torque of the isolated propeller are 5.2% lower and 9.1% higher respectively than the MRF-VLM results of VSPAERO. Although the MRF-VLM method has limitations of insufficient viscosity effect simulation, the design of the propeller involves multiple parameters, and the iteration can be completed within 5 min in the conceptual design stage. Compared with the calculation time of 1 h for a single working condition of MFR-RANS, the calculation efficiency has increased by more than 10 times.

STAR-CCM+ conducted MRF-RANS numerical simulations for the three configurations, including baseline WB, optimized WB and optimized WBP. Among them, the leading edge, trailing edge of the wing and the propeller area all adopted refinement processing. The total number of grids WB is 10.5 million, and the grids reached 14.58 million with adding the propeller rotation domain. The surface grid divisions of the WB and WBP are shown in Figures 18a,b, and the rotation domain and far-field surface grid are shown in Figures 18b,c.

FIGURE 18

In Figure 19, the STAR-CCM+ computational results reveal the significant impact of propeller slipstream on the aerodynamic characteristics of the wing. Compared to the baseline WB, optimized WB exhibits an expanded leading-edge negative pressure extending chordwise region. When the propeller slipstream is further introduced in the optimized WBP configuration, the leading-edge negative pressure zone expands even more. This trend of gradual negative pressure region enlargement is consistent with the evolutionary patterns of pressure difference contour plots observed during two optimization stages: WB optimization (Figure 12a,b) and WBP coupling optimization (Figure 15) based on the VLM.

FIGURE 19

The quantitative data in Table 7 shows that the C_l of the WBP configuration calculated by STAR-CCM+ increases by 15.6%, and the L/D increases by 18.4%, which is consistent with the trend of slipstream enhancement predicted by VSPAERO. Although the C_l errors between the two solvers is 9.7%, both effectively capture the slipstream enhancement on aerodynamic performance. Aerodynamic curves in Figure 20 further demonstrate that while numerical discrepancies exist in C_l and L/D predictions across different angles of attack, the variation trends in aerodynamic characteristics show strong agreement between STAR-CCM+ and VSPAERO. The calculation time of this high-fidelity RANS method for a single working condition of the WB model and the WBP model is 1 h and 1 h 20 min respectively. Compared with this method, the calculation efficiency of the VLM method has increased by 80–160 times.

FIGURE 20

Conclusion

In conclusion, the proposed method not only plays a significant role in the current eVTOL aerodynamic design but also lays a foundation for future applications in urban air traffic, where iterative exploration and physical feasibility are of crucial importance.

References

• 1

  Patterson MD Antcliff KR Kohlman LW . A Proposed Approach to Studying Urban Air Mobility Missions Including an Initial Exploration of Mission Requirements. In: Annual Forum and Technology Display (2018).

  • Google Scholar

• 2

  Roy S Maheshwari A Crossley WA DeLaurentis DA . Future Regional Air Mobility Analysis Using Conventional, Electric, and Autonomous Vehicles. J Air Transportation (2021) 29(3):113–26. 10.2514/1.D0235

  • CrossRef
  • Google Scholar

• 3

  Daskilewicz M German B Warren M Garrow LA Boddupalli SS Douthat TH . Progress in Vertiport Placement and Estimating Aircraft Range Requirements for EVTOL Daily Commuting. 2018 Aviation Technology, Integration, and Operations Conference. Atlanta, Georgia: AIAA2018–2884. (2018). 10.2514/6.2018-2884

  • CrossRef
  • Google Scholar

• 4

  Brown A Harris WL . Vehicle Design and Optimization Model for Urban Air Mobility. J Aircraft (2020) 57(6):1003–13. 10.2514/1.C035756

  • CrossRef
  • Google Scholar

• 5

  Brelje BJ Martins JRRA . Electric, Hybrid, and Turboelectric fixed-wing Aircraft: A Review of Concepts, Models, and Design Approaches. Prog Aerospace Sci (2019) 104:1–19. 10.1016/j.paerosci.2018.06.004

  • CrossRef
  • Google Scholar

• 6

  Ronzheimer A . Aircraft Geometry Parameterization with high-end CAD-Software for Design Optimization. In: Proceedings ECCOMAS (2012).

  • Google Scholar

• 7

  Bhise TA Course MD . Conceptual Design and Parametric Modeling of the Aircraft Wing (2020). 10.13140/RG.2.2.29170.94400

  • CrossRef
  • Google Scholar

• 8

  Lyu Z Martins JRRA . Aerodynamic Design Optimization Studies of a blended-wing-body Aircraft. J Aircraft (2014) 51(5):1604–17. 10.2514/1.C032491

  • CrossRef
  • Google Scholar

• 9

  Li L Bai J Guo T . Aerodynamic Optimization Design of the Supersonic Aircraft Based on Discrete Adjoint Method. J Northwest Polytechnical Univ (2017) 35(5):843–9. 10.3390/aerospace10050420

  • CrossRef
  • Google Scholar

• 10

  Bryant AH Walters DK . Developing a Small-Scale Propeller Thrust Model Using Experimentation and CFD. AIAA Scitech 2021 Forum. AIAA2021–1307. Virtual Event. (2021) 10.2514/6.2021-1307

  • CrossRef
  • Google Scholar

• 11

  Liu X Zhao D Oo NL . Comparison Studies on Aerodynamic Performances of a Rotating Propeller for small-size Uavs. Aerospace Sci Technology (2023) 133:108148. 10.1016/j.ast.2023.108148

  • CrossRef
  • Google Scholar

• 12

  Sinnige T van Arnhem N Stokkermans TCA Eitelberg G Veldhuis LLM . Wingtip-Mounted Propellers: Aerodynamic Analysis of Interaction Effects and Comparison with Conventional Layout. J Aircraft (2019) 56(1):295–312. 10.2514/1.c034978

  • CrossRef
  • Google Scholar

• 13

  Li P Zhao QJ . CFD Calculations on the Interaction Flow Field and Aerodynamic Force of tiltrotor/wing in Hover. Acta Aeronautica et Astronautica Sinica (2014) 5(2):361–71. 10.7527/S1000-6893.2013.0220

  • CrossRef
  • Google Scholar

• 14

  Xu J Yu J Lu X Long Z Xu Y Sun H . Aerodynamic Performance and Numerical Analysis of the Coaxial Contra-rotating Propeller Lift System in Evtol Vehicles. Mathematics (2024) 12(7):1056. 10.3390/math12071056

  • CrossRef
  • Google Scholar

• 15

  Dang TQ . Simulations of propeller/airframe Interference Effects Using an Euler Correction Method. J Aircraft (1989) 26(11):994–1001. 10.2514/3.45872

  • CrossRef
  • Google Scholar

• 16

  Thom A Duraisamy K . Computational Investigation of Unsteadiness in Propeller Wake–Wing Interactions. J Aircraft (2013) 50(3):985–8. 10.2514/1.C031866

  • CrossRef
  • Google Scholar

• 17

  Whitfield DL Jameson A . Euler Equation Simulation of propeller-wing Interaction in Transonic Flow. J Aircraft (1984) 21(11):835–9. 10.2514/3.45052

  • CrossRef
  • Google Scholar

• 18

  Benaouali A Kachel S. A Surrogate-based Integrated Framework for the Aerodynamic Design Optimization of a Subsonic Wing Planform Shape. Proc Inst Mech Eng G: J Aerospace Eng (2018) 232(5):872–83. 10.1177/0954410017699007

  • CrossRef
  • Google Scholar

• 19

  Rao KS Abhilasha AN Das A Sivapragasam M . Aerodynamic Shape Optimization at Low Reynolds Number Using multi-level Hierarchical Kriging Models. Optimization Eng (2025) 26(1):479–505. 10.1007/s11081-024-09915-2

  • CrossRef
  • Google Scholar

• 20

  Piotrowski MG Zingg DW . Investigation of a Smooth Local Correlation-Based Transition Model in a Discrete-Adjoint Aerodynamic Shape Optimization Algorithm. AIAA Scitech 2022 Forum. San Diego, CA & Virtual: AIAA2022–1865. (2022) 10.2514/6.2022-1865

  • CrossRef
  • Google Scholar

• 21

  Zhang W Tan M Liu YF Nie YB Luan Y . Wing Design and Adaptive Optimization of a Compound Drone. Phys Gases (2024) 9(2):54–65. (in Chinese). 10.19527/j.cnki.2096-1642.1091

  • CrossRef
  • Google Scholar

• 22

  Abdul-Kaiyoom MAS Yildirim A Martins JRRA . Coupled Aeropropulsive Design Optimization of an Over-wing Nacelle Configuration. J Aircraft (2025) 62(1):94–116. 10.2514/1.C037678

  • CrossRef
  • Google Scholar

• 23

  Chauhan SS Martins JRRA . Rans-Based Aerodynamic Shape Optimization of a Wing Considering propeller–wing Interaction. J Aircraft (2021) 58(3):497–513. 10.2514/1.C035991

  • CrossRef
  • Google Scholar

• 24

  Pedreiro LN . Estudo E Otimização De Uma Asa Sob Efeito De Hélice Na Configuração Tractor Para Redução De Arrasto (2017).

  • Google Scholar

• 25

  Alba C Elham A German BJ Veldhuis LL . A Surrogate-based multi-disciplinary Design Optimization Framework Modeling wing–propeller Interaction. Aerospace Sci Technology (2018) 78:721–33. 10.1016/j.ast.2018.05.002

  • CrossRef
  • Google Scholar

• 26

  Paek SK Chae S Kang HJ . Development of Mission Analysis Tool for Evtol Aircrafts of Lift-Tilt Concept. J Korean Soc Aeronaut & Space Sci (2021) 49(10):859–70. 10.5139/JKSAS.2021.49.10.859

  • CrossRef
  • Google Scholar

• 27

  Sun ZY Wang Q Qiao W Fu C Zheng YF Song JY . Research on the Analysis Method of Distributed propeller-wing Design of Electric Aircraft. Aeronaut Sci & Technology (2024) 35(1):53–64 (in Chinese). 10.19452/j.issn1007-5453.2024.01.006

  • CrossRef
  • Google Scholar

• 28

  Yu HL Lei T Zhang XY Zhang XB . A Fast Multidisciplinary Parameter Estimation Method for Distributed Electric Propulsion Aircraft. Aeronaut Sci Technolog (2024) 35(1):65–74 (in Chinese). 10.19452/j.issn1007-5453.2024.01.007

  • CrossRef
  • Google Scholar

• 29

  Sheridan CN Pham DD Whiteside S . Evaluation of Vspaero Analysis Capabilities for Conceptual Design of Aircraft with Propeller-Blown Wings. (2021). AIAA 2021–2510. Virtual Event. 10.2514/6.2021-2510

  • CrossRef
  • Google Scholar

• 30

  Wu J Gao F Li S Yang F . Conceptual Design and Optimization of Distributed Electric Propulsion General Aviation Aircraft. Aerospace (2023) 10(5): 387. 10.3390/aerospace10050387

  • CrossRef
  • Google Scholar

• 31

  Anderson JD . Fundamentals of Aerodynamics(Seventh Edition). New York: McGraw Hill Higher Education (2024).

  • Google Scholar

• 32

  Katz J Plotkin A . Low-Speed Aerodynamics. Cambridge: Cambridge University Press (2001).

  • Google Scholar

• 33

  Drela M . Flight Vehicle Aerodynamics. MIT press (2014).

  • Google Scholar

• 34

  Cheng ZY Yang YX Zhang XC Yu LF Ye B . Rapid Evaluation Method for Aerodynamic Characteristics of Distributed Electric Propulsion Aircraft Concept Scheme. J Beijing Univ Aeronautics Astronautics (2023) 49(11):3047–58. (in Chinese). 10.13700/j.bh.1001-5965.2021.0771

  • CrossRef
  • Google Scholar

• 35

  Miranda LR Elliot RD Baker WM . A Generalized Vortex Lattice Method for Subsonic and Supersonic Flow Applications. (1977).

  • Google Scholar

• 36

  Resende GJ Malatesta V Savio MC Castro BM . Wing’s Aerodynamic Characteristics due to Distributed Propulsion over the Wingspan. J Braz Soc Mech Sci Eng (2023) 45(9):495. 10.1007/s40430-023-04413-5

  • CrossRef
  • Google Scholar

• 37

  Gudmundsson S . General Aviation Aircraft Design: Applied Methods and Procedures. Butterworth-Heinemann (2013).

  • Google Scholar

• 38

  Conway JT . Analytical Solutions for the Actuator Disk with Variable Radial Distribution of Load. J Fluid Mech (1995) 297:327–55. 10.1017/s0022112095003120

  • CrossRef
  • Google Scholar

• 39

  Johnson W . Helicopter Theory. Princeton: Princeton University Press (1980).

  • Google Scholar

• 40

  McAlister KW . NACA 0015 Wing Pressure and Trailing Vortex Measurements. Technical Paper of National Aeronautics and Space Administration, Office of Management, Scientific and Technical Information Program. NASA-TP-3151. (1991).

  • Google Scholar

• 41

  Veldhuis LLM . Propeller Wing Aerodynamic Interference. Delft: Delft University of Technology (2005).

  • Google Scholar

• 42

  Jin R Chen W Sudjianto A . An Efficient Algorithm for Constructing Optimal Design of Computer Experiments. J Stat Plann inference (2005) 134(1):268–87. 10.1016/j.jspi.2004.02.014

  • CrossRef
  • Google Scholar

• 43

  Sacks J Welch WJ Michell TJ Wynn HP . Design and Analysis of Computer Experiments. Stat Sci (1989) 4(4):409–35. 10.1214/ss/1177012413

  • CrossRef
  • Google Scholar

• 44

  Zhu XQ . Research on Flight Vehicle Shape Optimization Design Algorithms and multi-output Surrogate Model. Xi’an: Northwestern Polytechnical University (2019).

  • Google Scholar

• 45

  Zhang J Sanderson AC . JADE: Adaptive Differential Evolution with Optional External Archive. IEEE Trans Evol Comput (2009) 13(5):945–58. 10.1109/TEVC.2009.2014613

  • CrossRef
  • Google Scholar

• 46

  Guo SM Yang CC Hsu PH Tsai JSH . Improving Differential Evolution with a successful-parent-selecting Framework. IEEE Trans Evol Comput (2015) 19(5):717–30. 10.1109/TEVC.2014.2375933

  • CrossRef
  • Google Scholar

• 47

  Guo SM Yang CC . Enhancing Differential Evolution Utilizing Eigenvector-based Crossover Operator. IEEE Trans Evol Comput (2014) 19(1):31–49. 10.1109/TEVC.2013.2297160

  • CrossRef
  • Google Scholar

• 48

  Wang Y Li HX Huang T Li L . Differential Evolution Based on Covariance Matrix Learning and Bimodal Distribution Parameter Setting. Appl Soft Comput (2014) 18:232–47. 10.1016/j.asoc.2014.01.038

  • CrossRef
  • Google Scholar

• 49

  Jones DR . A Taxonomy of Global Optimization Methods Based on Response Surfaces. J Glob optimization (2001) 21:345–83. 10.1023/a:1012771025575

  • CrossRef
  • Google Scholar

• 50

  Kim DH Lee JH Hwang HY . Aerodynamic Analysis, Required Power and Weight Estimation of a Compound (Tilt Rotor+ Lift+ Cruise) Type Evtol for Urban Air Mobility Using Reverse Engineering Techniques. J Adv Navigation Technology (2021) 25(1):17–28. 10.12673/JANT.2021.25.1.17.c

  • CrossRef
  • Google Scholar