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: V → P = Γ 4 π · ∫ a b d l → × r → − r → ′ r p 3 (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 β = 1 − M ∞ 2 is introduced to account for compressibility effects [33]. Induced velocities are computed using the Biot-Savart law, expressed as Equation 2: V → P = − β 2 π K · ∫ a b Γ → × r → − r → ′ r β 3 d l (2) where Γ → represents the strength of the vortex ring, and K is the factor considering the influence of compressibility on the calculation: if M ∞ <   1 , then K = 2; otherwise, K = 1, r β is defined by Equation 3 r β 2 = x − x ′ 2 + β 2 y − y ′ 2 + z − z ′ 2 (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].

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: V z r , z = 2 V z r , 0 + v i · − α + z R arcsin k 0 z ≥ 0 (4) V z r , z = v i · α + z R arcsin k 0 z < 0 (5)

The radial velocity is given by Equation 6, k 0 is defined by Equation 7, and the pressure jump across the disk surface is given by Equation 8. V r r , z = v i z 2 r 1 α − α − v i · r 2 R arcsin k 0 (6) k 0 = 2 R z 2 + R + r 2 1 / 2 + z 2 + R − r 2 1 / 2 (7) Δ P r = 2 ρ ∞ V ∞ V z r , 0 (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, V z r , 0 = v i R R 2 − r 2 1 / 2 (9) v i = − V ∞ 2 + V ∞ 2 2 + T 2 ρ ∞ A (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. α = R a 2 − r 2 − z 2 2 + 4 R a 2 z 2 1 / 2 + R a 2 − r 2 − z 2 2 R a 2 1 / 2 (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 V ∞ is defined by Equation 13. V t r = V ∞ + v o v o ω r ω r 2 + V ∞ + v o 2 + 2 v o C d C l (12) v o 2 = v i 2 / 1 + C T ∼ ⁡ ln C T ∼ + C T ∼ / 2 (13)

The C _l /C _d is modeled by combining the blade hover thrust coefficient ( C T ∼ = T / ρ ∞ A ω R 2 ) and the hover power coefficient ( C P ∼ = P / ρ ∞ A ω R 3 ), leading to Equations 14, 15. σ C l = 6 C T ∼ (14) C P ∼ = σ C d 8 + κ C T ∼ 3 / 2 2 (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, C T ∼ and C P ∼ 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.

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 ( V ∞ = 44.6 m/s), Reynolds number Re = 1.5 × 10⁶, and air density ρ = 1.225 kg/m³.

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.

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.

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 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.

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.

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.

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.

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.

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.

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.

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.

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. minimize − L / D w . r . t . x ∈ R 23 subject to C l w i n g ≥ C l w i n g , 0 0.9 S w i n g , 0 ≤ S w i n g ≤ 1.1 S w i n g , 0 (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.

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.

With the optimized drag coefficient Cd = 0.0273 and drag Equation 17, the total cruise drag F is calculated as 935 N. F = 1 2 ρ v 2 C d S wing (17)

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. T t o t a l = F η t o t a l (18)

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 _P Equation 19 is inversely deduced by the efficiency Equation. C P ≥ C T J η max η max = 0.9 (19)

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 minimize − η * L / D w . r . t . x ∈ R 7 subject to C T ∈ 0.02 , 0.06 C P ∈ 0.04 , 0.08 η ∈ 0.75 , 0.9 (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.

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. V s l i p s t r e a m = V ∞ 1 + C T 2 (21)

The increase in diameter reduces both C _T and slipstream velocity V s l i p s t r e a m . 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.

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.

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.

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.

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.

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.

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.