© Md. Hakim et al., Published by EDP Sciences, 2026

1 Introduction

Mini-unmanned aerial vehicles have gained increasing interest for a wide range of civil missions, including aerial imaging, communications, agricultural operations, and light-delivery tasks. In comparison with piloted aircraft, mini-UAVs offer greater layout flexibility because they are not constrained by cockpit integration or direct pilot accommodation. At the same time, their mission dependent requirements and reduced scale make the direct application of conventional aircraft conceptual design practices less straightforward [1,2]. This context has encouraged the exploration of unconventional airframe architectures capable of achieving compactness, aerodynamic efficiency, and functional integration. Among these, BWB concept has emerged as a promising candidate for mini-UAV applications [3–5].

BWB configurations combine lifting and volume-carrying functions within a highly integrated airframe, which can provide aerodynamic and structural advantages over conventional layouts. These benefits include improved lift-to-drag characteristics, reduced wetted area, and favorable internal volume distribution [6]. However, when the configuration is tailless, these advantages are accompanied by a major design challenge where the aircraft must achieve satisfactory static stability without relying on conventional horizontal and vertical tail surfaces. As a result, the restoring moments required for stable flight must be generated primarily through wing geometry, airfoil characteristics, twist distribution, and center-of-gravity placement [7]. For tailless BWB mini-UAVs, static stability is therefore not a secondary design check, but a main constraint that must be addressed from the earliest conceptual-design stage.

This challenge is particularly important because the longitudinal and lateral-directional stability characteristics of tailless aircraft depend on several coupled geometric parameters. Variations in sweep, chord distribution, twist, dihedral, and planform break location can significantly alter the aerodynamic derivatives governing pitch, roll, and yaw behavior. A meaningful exploration of these effects requires the evaluation of many candidate geometries across a multidimensional design space. Yet this rapidly becomes expensive when based solely on high-fidelity numerical tools. Although Computational Fluid Dynamics (CFD) can provide detailed aerodynamic information, its computational cost makes large parametric studies less practical during early design stages, where rapid screening and iterative refinement are required [8,9].

To address this limitation, surrogate modeling has become an attractive tool in aerospace design. Surrogate models, also referred to as metamodels, approximate the response of a numerical model by learning the relationship between design inputs and output quantities from a limited set of simulations [10]. Techniques such as polynomial regression, Kriging, and Gaussian Process Regression have been widely used to support optimization, sensitivity analysis, uncertainty quantification, and design-space exploration while reducing computational cost [11–13]. In aerodynamic applications, these methods are especially valuable when repeated evaluations are needed, as in conceptual design studies involving many geometric configurations.

Several previous works have demonstrated the usefulness of surrogate-based methods in aerodynamic design. Jeong et al. [14] applied a Kriging-based functional analysis of variance to two-dimensional airfoil optimization and identified the most influential design variables. Norazila and Masahiro [15] developed a surrogate framework for aerodynamic coefficients and derivatives intended for flight simulation. Zhang et al. [9] proposed a collaborative surrogate-modeling method using multi-fidelity aerodynamic data, while Andrés-Pérez and Paulete-Periáñez [16] assessed regression-based surrogates as alternatives to more expensive CFD calculations for aerodynamic-coefficient estimation. Taken together, these studies confirm that surrogate techniques can provide substantial computational savings while preserving adequate predictive capability for early-stage aerodynamic investigations.

From the stability point of view, previous studies have emphasized the strong sensitivity of tailless configurations to geometric design choices. Shams et al. [17] showed that parameters such as sweep, dihedral, and aspect ratio significantly affect the aerodynamic and stability characteristics of flying-wing micro air vehicles operating at low Reynolds numbers. These results are particularly relevant because they show that, for tailless configurations, relatively small geometric variations can significantly modify the stability behavior and must therefore be considered carefully during early-stage design.

Surrogate-assisted design has also been applied to unconventional UAV aerodynamic optimization. Nikolaou et al. [18] presented a surrogate-modeling framework for winglet design and aerodynamic-performance optimization of UAVs, showing the practical value of data-driven methods for rapid exploration of geometric design alternatives. Although this study is not specifically focused on BWB stability, it further demonstrates that surrogate-based methodologies can efficiently support conceptual and preliminary aerodynamic design of unconventional UAV configurations.

Blended-wing-body configurations have received increasing attention in recent years from both stability and conceptual-design perspectives. Bainsla et al. [19] investigated the design, control-surface optimization, and stability analysis of a BWB unmanned aerial vehicle, highlighting the importance of stability considerations in such highly integrated aircraft. Wauters [20] later examined design optimization-under-uncertainty for a forward-swept BWB UAV using a surrogate-assisted framework, showing how computationally efficient methodologies can be used to satisfy aerodynamic and stability requirements simultaneously. More recently, Lyu et al. [21] proposed a nonlinear multi-fidelity aerodynamic surrogate model for blended-wing-body aircraft conceptual design optimization, further illustrating the growing relevance of surrogate-based approaches in BWB design-space exploration.

Despite this progress, the literature still offers limited work on surrogate modeling specifically dedicated to the static-stability assessment of tailless BWB mini-UAVs. Existing studies have generally focused on broader aerodynamic design, control-surface optimization, uncertainty-aware configuration design, or overall performance improvement, while stability-oriented works often do not construct a unified surrogate framework for the direct prediction and screening of the key longitudinal and lateral-directional static-stability derivatives across a parameterized tailless BWB design space. This is therefore a meaningful gap at the conceptual-design level, where designers need a computationally efficient method capable of rapidly identifying geometries that satisfy the main static-stability requirements before moving to higher-fidelity CFD or experimental validation.

The present work addresses this need by proposing an integrated framework for the static-stability design of tailless BWB mini-UAVs based on XFLR5 aerodynamic analysis, Latin Hypercube Sampling (LHS), and Gaussian Process Regression. In the proposed approach, XFLR5 is used to generate the aerodynamic database associated with a parameterized family of BWB geometries, LHS is employed to explore the design space efficiently, and GPR is used to construct surrogate models for the principal static-stability metrics. This combination enables rapid prediction of the aerodynamic derivatives governing longitudinal and lateral-directional stability while maintaining a computational cost compatible with conceptual design studies.

The main contribution of this study is therefore not surrogate modeling in isolation, but the development of a unified surrogate-based framework for static-stability screening and interpretation in tailless BWB mini-UAVs. Beyond identifying feasible statically stable configurations, the method is used to clarify the relative influence of selected geometric variables on the principal stability derivatives, thereby providing practical guidance for early-stage design decisions.

The remainder of this paper is organized as follows. Section 2 presents the flight-mechanics background and the static-stability criteria used in this work, together with the geometric parameterization of the BWB configuration. Section 3 justifies the aerodynamic modeling strategy and the choice of XFLR5 for conceptual-design analysis. Section 4 describes the surrogate-modeling framework, including Gaussian Process interpolation, Latin Hypercube Sampling, and cross-validation. Section 5 presents the model-selection results and the sensitivity analysis of the stability derivatives. Section 6 discusses the feasible BWB configurations identified through the proposed framework. Finally, Section 7 summarizes the main conclusions and outlines directions for future work.

2 Flight mechanics background and static stability formulation

The design of any aircraft, including the BWB UAV, fundamentally relies on flight mechanics modeling and analysis. These analyses form the foundation for evaluating performance, controllability, and stability characteristics throughout the conceptual and preliminary design stages. The classical theory of flight dynamics provides a rigorous framework for describing the motion of an aircraft under the action of aerodynamic, gravitational, and control forces. Comprehensive derivations of the governing equations can be found in standard references such [7,22–25], which together establish the theoretical background for aircraft stability and control.

According to Roskam [24], static stability refers to the inherent tendency of an aircraft to generate restoring forces and moments that oppose any small disturbance from an equilibrium or trimmed condition. These perturbations may arise from atmospheric turbulence, control inputs, or other transient external effects. When an aircraft is statically stable, a displacement from its trim state induces aerodynamic forces and moments that naturally drive it back toward that original state without continuous pilot or autopilot intervention. This concept of initial tendency distinguishes static stability from dynamic stability, which concerns the time-dependent response following a disturbance.

In the remainder of this section, we summarize the derivations of relevant stability characteristics governing the aircraft’s ability to restore the original trimmed conditions. Figure 1 shows the body-axis system (x_b,y_b,z_b); with origin at center of mass of the vehicle. The (X, Y, Z) and (L, M, N) are respectively, aerodynamic forces and moments. (u, v, w) and (p, q, r) are respectively linear and rotational velocities. These later define rotational motions in roll (p), pitch (q) and yaw (r).

Thrust loads are not invoked in this study for two main reasons; first is that, no engine failure is considered, and second, their derivatives are relatively small, so only the aerodynamics are prevalent [7]. They still, of course, of major importance in trim conditions, controls and maneuverings.

As generally conceived, x - z is a plane of symmetry, so symmetric or longitudinal motion; such in rectilinear steady state flight, is frequently separated from the coupled asymmetric lateral/directional motion; as in level turn or under crosswind [7]. Wind or stability axis is another option preferred for aerodynamic load computations [25].

Fig. 1 Body axis system and motion variables.

2.1 Longitudinal static stability

Among the various aspects of static stability, longitudinal static stability plays a key role in an aircraft’s ability to maintain a steady level flight path. It reflects how the vehicle responds to changes in its angle of attack (AoA) and directly governs its pitching behavior. Figure 2 shows the main aerodynamic forces and moments acting on an aircraft during steady symmetric flight, including drag (D), lift (L), and weight (W), as well as the AoA (α), and the moment about the center of gravity (CG).

Longitudinal static stability ensures that when a disturbance causes a small increase in angle of attack, the aircraft naturally produces a nose-down pitching moment that tends to restore the original equilibrium. Conversely, when the angle of attack decreases, a nose-up moment is generated to reestablish trim. This self-correcting tendency allows the aircraft to maintain a stable trajectory in the presence of small perturbations, ensuring predictable and controllable flight characteristics without requiring continuous control input.

In tailless or blended-wing configurations, longitudinal stability becomes even more critical, as the absence of a horizontal tailplane limits the designer’s ability to balance pitching moments through conventional tail lift. Instead, stability must be achieved through careful selection of the aerodynamic center location, airfoil pitching moment, planform twist distribution, and center of gravity placement.

This type of stability governs the aircraft’s response to changes in the angle of attack, which directly affects the pitching moment. In the case of longitudinal stability, an increase in the angle of attack should induce a negative (nose-down) pitching moment, reducing the angle of attack and bringing the aircraft back to equilibrium. The converse also holds true for a decrease in AoA. Proper longitudinal stability is crucial for ensuring that the aircraft maintains a stable trajectory and is able to recover from small perturbations in flight, ensuring predictable and controllable behavior.

For a conventional wing-body system, the governing expression of longitudinal static stability arises from the pitching-moment equilibrium condition, as detailed by Nelson [7]. The pitching moment coefficient about the center of gravity, C_mcg, can be expressed as:

$$C_{m_{cg}}=C_{m_{ac}}+C_{L_0}\left(\frac{x_{cg}}{\stackrel{‐}{c}}-\frac{x_{ac}}{\stackrel{‐}{c}}\right)+C_{L_{\alpha }}\left(\frac{x_{cg}}{\stackrel{‐}{c}}-\frac{x_{ac}}{\stackrel{‐}{c}}\right)\alpha$$(1)

Where $C_{m_{cg}},C_{m_{ac}},C_{L_0},C_{L_{\alpha }},x_{cg},x_{ac},\alpha$ and ^--c̄ are respectively the pitching moment coefficient at cg, the pitching moment coefficient at aerodynamic center (ac), lift coefficient at zero angle of attack, the lift curve slope, center of gravity position, aerodynamic center position (neutral point), angle of attack and the mean chord.

From equation (1), the classical longitudinal static stability conditions for a wing-only or tailless configuration are derived. A BWB aircraft is statically stable if the aerodynamic center lies aft of the center of gravity and the aircraft is trimmed at a positive angle of attack. These requirements lead to the following two stability criteria:

$$\frac{\partial C_m}{\partial \alpha }=C_{m_{\alpha }}=C_{L_{\alpha }}\left(\frac{x_{cg}}{\stackrel{‐}{c}}-\frac{x_{ac}}{\stackrel{‐}{c}}\right)<0$$(2)

$$C_{m_0}=C_{m_{ac}}+C_{L_0}\left(\frac{x_{cg}}{\stackrel{‐}{c}}-\frac{x_{ac}}{\stackrel{‐}{c}}\right)\ge 0.$$(3)

Equation (2) indicates that the slope of the pitching-moment curve with respect to the angle of attack must be negative for stability, implying that an increase in α produces a restoring (nose-down) moment. Equation (3) expresses the trim condition, requiring that the aircraft’s pitching moment at zero lift be positive or zero to permit a stable equilibrium at a positive angle of attack.

For more detailed evaluation, the pitching moment coefficient about the aerodynamic center, C_mac, can be determined from the spanwise integration of local airfoil characteristics and planform geometry as:

$$C_{m_{ac}}=\frac{1}{S\stackrel{‐}{c}}\underset{-\frac{b}{2}}{\overset{\frac{b}{2}}{{\displaystyle \int }}}\left(c_{m_{ac}}\left(y\right)c{\left(y\right)}^2\right)dy+\pi \underset{-\frac{b}{2}}{\overset{\frac{b}{2}}{{\displaystyle \int }}}\left({\alpha }_{0_L}+ϵ\left(y\right)-{\alpha }_0\left(y\right)\right)c\left(y\right)x_1\left(y\right)dy$$(4)

Where S is the surface, c_mac is the airfoil pitching moment, c is the chord, α_0L is angle of attack for zero planform lift, ε is the planform twist angle, α₀ is section airfoil zero lift angle of attack, x₁ the x displacement between the quarter chord points in the mean aerodynamic and the actual chord at a given spanwise location y on the wing and b is the wingspan.

Fig. 2 Aerodynamic forces and pitching moment about center of gravity.

2.2 Lateral-directional static stability

In addition to longitudinal stability, an aircraft must also maintain strong lateral–directional static stability to ensure safe and coordinated flight. These two coupled modes of stability determine how well the aircraft can hold its attitude and heading when faced with asymmetric disturbances such as gusts, crosswinds, or uneven lift distribution.

Lateral stability describes the aircraft’s tendency to resist rolling and to return to wings level flight after a disturbance in bank angle. Directional stability, on the other hand, relates to the aircraft’s ability to resist yawing and maintain its intended heading (see Fig. 3). These characteristics allow the aircraft to recover smoothly from minor deviations and continue flying in a coordinated manner without requiring excessive pilot input or automated control corrections.

A well-tuned lateral–directional stability ensures that disturbances don’t develop into unwanted motions like spiral divergence or Dutch roll. Achieving this balance depends on a several geometric and aerodynamic factors, mainly wing sweep, dihedral angle, vertical tail size, and the lateral distribution of lift and drag. In conventional aircraft, these stabilizing effects are provided by the horizontal and vertical tails. However, in tailless configurations, produce the same stabilizing moments demands precise aerodynamic shaping and twist distribution.

Applying the principle of static stability to rotation about the yaw axis implies that a stable airplane must demonstrate weathercock stability. When the aircraft experiences a small sideslip angle, β, relative to its flight path, the resulting aerodynamic forces should produce a yawing moment that acts to restore it to symmetric flight. For this restoring tendency to occur, the derivative of the yawing-moment coefficient with respect to sideslip angle must be positive [7]:

$$C_{N_{\text{β}}}=\frac{d{C}_N}{d\beta }>0.$$(5)

A positive C_Nβ indicates that when the aircraft yaws relative to the oncoming airflow, the aerodynamic forces generate a corrective moment that aligns the nose back into the relative wind, thus maintaining directional stability.

Lateral static stability (roll stability), on the other hand, is associated with the restoring rolling moment produced when the aircraft is displaced from a wings-level attitude. When a sideslip occurs, the differential lift generated across the wings induces a rolling moment that can either amplify or counteract the bank angle. A laterally stable configuration produces a restoring (opposite) roll moment that tends to relevel the wings. Mathematically, the criterion for static roll stability can be expressed as [7]:

$$C_{l_{\beta }}=\frac{d{C}_l}{d\beta }<0.$$(6)

A negative C_lβ ensures that the aircraft will roll in the opposite direction of the sideslip, thus promoting automatic recovery to a symmetric attitude.

Several studies [7, 22–24] have demonstrated that the rolling moment generated during a sideslip is influenced by factors such as wing dihedral, sweep angle, wing placement on the fuselage, and the presence of a vertical tail. In the case of a BWB, the latter two stabilizing features are absent.

Fig. 3 Aerodynamic forces and pitching moment about center of gravity.

2.3 Design variables for static stability analysis of BWB UAV

The inherent static stability and overall aerodynamic behavior of a BWB UAV are strongly governed by its geometric configuration. In tailless aircraft, where no horizontal or vertical tail is available to generate conventional restoring moments, stability must be achieved primarily through the shaping of the wing itself. As a result, relatively small variations in planform geometry, twist distribution, or dihedral can significantly affect the aerodynamic derivatives governing pitch, roll, and yaw behavior. As summarized in Table 1, the key static-stability derivatives considered in this work are influenced by geometric parameters such as sweep, twist, chord distribution, and dihedral.

The seven design variables retained in this study were selected based on a review of the literature on tailless, flying-wing, and blended-wing-body aircraft, where sweep, twist, chord distribution, and dihedral are consistently identified as the main geometric parameters affecting aerodynamic and static-stability characteristics. The selected variables were therefore not chosen arbitrarily, but were defined so as to capture the dominant geometric effects reported in previous studies while keeping the parameterization sufficiently compact for surrogate-model construction at the conceptual-design stage.

More specifically, the root chord, kink chord, kink offset, and kink sweep were selected to describe the main features of the wing planform and its spanwise loading distribution. The root-tip twist difference was retained because twist directly affects the pitching-moment behavior and the longitudinal lift distribution, which are essential for trim and longitudinal static stability in tailless configurations. The inward and outward dihedral angles were included because they influence the rolling and yawing moments generated in sideslip and therefore play a direct role in lateral-directional stability. Together, these seven variables provide a physically meaningful yet manageable description of the geometric mechanisms that govern the stability behavior of the considered BWB UAV.

Other potentially relevant parameters were kept fixed in the present work either because they were prescribed by the baseline configuration, regarded as secondary with respect to the objectives of the study, or excluded to avoid an excessively large design space that would reduce the interpretability and efficiency of the surrogate-modeling framework. In particular, the airfoil was selected beforehand based on the criteria and analysis presented in [26], with emphasis on a low negative pitching-moment coefficient about the aerodynamic center, a high maximum angle of attack, and low drag characteristics. These criteria are particularly appropriate for tailless and BWB applications, where longitudinal stability must be achieved without the benefit of a horizontal tail. Since the purpose of the present work is not to optimize the airfoil itself, but rather to investigate the influence of global geometric parameters on static stability, the airfoil is treated as fixed throughout the analysis.

The baseline configuration shown in Figure 4 corresponds to a pure flying-wing architecture with no conventional fuselage or tail assemblies. In this all-wing concept, aerodynamic and stability requirements must be satisfied solely through the geometry of the lifting surface and its twist distribution. For this reason, the adopted baseline is representative of a highly integrated UAV platform in which stability must be achieved inherently through wing geometry alone.

To systematically assess how geometric variations influence stability, the seven design variables were defined around this baseline configuration, as summarized in Table 2. Each variable was assigned a baseline value and a variation range chosen to capture its aerodynamic influence while preserving geometric feasibility and structural plausibility. These ranges provide sufficient coverage of the design space to generate representative aerodynamic data and to support the construction of reliable surrogate models for conceptual static-stability analysis.

Fig. 4 Baseline tailless BWB geometry and definition of the seven design variables used in the parametric study.

3 Choice of numerical model for stability analysis

In the preliminary design of unconventional configurations such as tailless BWB UAVs, the selection of an aerodynamic analysis method must balance predictive capability and computational efficiency. At this stage, the objective is not the final aerodynamic certification of a single geometry, but the rapid evaluation of many candidate configurations in order to understand how geometric variations affect the static-stability derivatives. The numerical tool must therefore be sufficiently accurate to capture the main aerodynamic trends, while remaining fast enough to support design-space exploration and surrogate-model construction.

High-fidelity CFD methods based on the Navier–Stokes equations, particularly Reynolds-Averaged Navier–Stokes (RANS) formulations, provide detailed information on viscous effects, boundary-layer development, and flow separation. For this reason, they are indispensable in later design stages and final aerodynamic validation. However, their use in conceptual design remains constrained by high computational cost and preprocessing effort. Generating suitable volume meshes, defining boundary conditions, and ensuring convergence and mesh independence require substantial time and computational resources, especially when a large number of geometric variants must be analyzed. As a result, CFD-based approaches are not the most practical option for rapid parametric studies at the early design stage [25].

For conceptual design purposes, potential-flow methods provide a more suitable compromise between cost and fidelity. Based on inviscid, incompressible, and irrotational flow assumptions, these methods can reproduce the global aerodynamic behavior of lifting configurations with much lower computational expense than CFD. Since they rely on surface discretization rather than full volume meshing, both preprocessing and solution time are greatly reduced. This makes them particularly attractive for design-of-experiments studies and surrogate-model generation, where a large number of simulations must be performed efficiently [25, 27].

In the present work, aerodynamic analyses were carried out using XFLR5 [28], an open-source tool widely used for low-speed aircraft design. XFLR5 combines XFOIL-based two-dimensional airfoil analysis with three-dimensional potential-flow solvers for complete configurations. The XFOIL module, originally developed by Drela [29], couples a panel method for the external inviscid flow with an integral boundary-layer formulation, allowing improved prediction of low-Reynolds-number aerodynamic behavior. For three-dimensional configurations, XFLR5 offers lifting-line, vortex-lattice, and panel methods. Among these, the three-dimensional panel method is the most suitable for the present BWB geometry, since it can represent thicker and more integrated lifting surfaces more effectively than simpler lifting-line or vortex-lattice approaches [28–30].

The main reason for selecting XFLR5 in this study is its suitability for conceptual static-stability analysis. It allows rapid computation of the aerodynamic forces, moments, and stability derivatives required to train the surrogate models, while maintaining a level of fidelity adequate for preliminary design assessment. In particular, it is well adapted to low-speed attached-flow conditions and to the comparative analysis of multiple candidate geometries. Compared with CFD solvers, XFLR5 requires orders of magnitude less computational effort, which is a key advantage when generating the aerodynamic database needed for Latin Hypercube Sampling and Gaussian Process Regression [30, 31].

It should be emphasized that the present stability analysis is restricted to linear static-stability derivatives evaluated in the vicinity of the reference operating condition. The framework therefore captures the initial restoring tendencies associated with small perturbations in angle of attack and sideslip, which is appropriate for conceptual static-stability assessment of tailless BWB configurations. However, as a potential-flow-based tool, XFLR5 does not fully resolve viscous effects, strong flow separation, stall onset, vortex-dominated interactions, or other nonlinear aerodynamic phenomena that may become important at higher angles of attack or in more complex low-Reynolds-number flow regimes. The present methodology should therefore be regarded as a conceptual-design and screening framework rather than a substitute for higher-fidelity CFD or experimental validation. Its role is to efficiently identify promising BWB geometries, reveal the main geometric trends governing static stability, and generate the aerodynamic database used for the surrogate-model construction described in the following section.

4 Gaussian process interpolation and surrogate construction

Following the aerodynamic analysis using XFLR5, the next stage of the study involves exploring the static-stability design space defined by the geometric variables introduced in Section 2.3. Since evaluating aerodynamic responses across a broad range of parameter combinations is time intensive using even potential-flow solvers such as XFLR5. This challenge becomes more pronounced when dozens or even hundreds of geometric configurations must be simulated to capture the relationships between input variables and the resulting stability derivatives. To reduce this computational load and facilitate systematic design exploration, surrogate modeling is introduced.

Surrogate models is a statistical approximation that emulates the input–output behavior of a computational model. It provides a means of estimating aerodynamic responses at new design points without the need to rerun the full aerodynamic simulation. In this work, the surrogate model is constructed to approximate the input–output relationships between the geometric parameters of the BWB configuration and the aerodynamic stability derivatives computed using XFLR5. Each design point represents a unique combination of design variables within the defined ranges (see Tab. 2), corresponding to one potential UAV configuration.

Among existing surrogate modeling approaches, Gaussian Process Regression has proven to be particularly powerful for aerospace design applications. Originally developed by Matheron [32] and later adapted for computer experiments by Sacks et al. [33], the Kriging method models the unknown system response as a realization of a stochastic process. This probabilistic formulation enables the model to interpolate the available data points exactly while also quantifying the uncertainty associated with its predictions. Gaussian Process Regression (GPR) is particularly effective when working with limited datasets or complex nonlinear relationships, as it offers an optimal balance between flexibility and predictive accuracy.

By integrating fast aerodynamic solvers such as XFLR5 with surrogate modeling techniques, it becomes possible to conduct extensive design studies, sensitivity analyses, and optimization tasks without incurring the high computational costs of full CFD simulations. This combined approach supports rapid iteration and informed decision-making early in the design process, prior to the more resource-intensive validation using high-fidelity solvers.

4.1 Gaussian process interpolation

In deterministic computer experiments, the simulation outputs are assumed to be noise-free, meaning that identical input conditions always yield identical results. Consequently, the surrogate model should interpolate the available data rather than approximate it with random error. Gaussian Process interpolation satisfies this requirement by representing the model response as a combination of a deterministic regression component and a stochastic process that captures the spatial correlation between nearby design points.

The observed computer outputs $Y=\left\{y\left(1\right),y\left(2\right),...,y\left(n\right)\right\}$ which represent evaluations of the unknown deterministic function y = f (x) at n sampled design point $x_i\left(x_{\text{i}1},\dots ,x_{ia}\right)$ where a is the number of design variables. The Gaussian Process model is expressed as

$$Y\left(x\right)={\displaystyle \sum }_{j=1}^p{f}_j\left(x\right){\text{β}}_j+Z\left(x\right)=F^T\left(x\right)\beta +Z\left(x\right)$$(7)

Where, β is the vector of regression coefficients to be estimated from the data, f_j(x) are the chosen basis function (or regressors), and Z(x) is a Gaussian stochastic process with zero mean and covariance that defines the spatial correlation of the function.

The covariance structure is given by:

$$Cov\left(Z\left(x\right),Z\left(x^{\text{'}}\right)\right)={\sigma }_z^{2}R\left(x,x\text{'}\right)$$(8)

Where ${\sigma }_z^2$ is the process variance and $R\left(\cdot ,\cdot \right)$ is a parametrized correlation function describing the dependency between two input points. A commonly adopted formulation is the Gaussian correlation function:

$$R\left(x,x\text{'}\right)={\displaystyle \prod }_{j=1}^a\exp \left(-{\theta }_j{\left|x_i-x_j^{\text{'}}\right|}^m\right)$$(9)

Where θ_j ≥ 0, and 0 < m ≤ 2 are hyperparameters controlling the correlation decay and the smoothness of the response. For highly differentiable aerodynamic functions, m = 2 is generally used.

The Best Linear Unbiased Predictor (BLUP) of the Gaussian Process is obtained by maximizing the likelihood function (l) of the observed data with respect to the model hyperparameters. The predictive mean at a test point (x^te) is given by [29,30]:

$$Y\left(x^{te}\right)=Y^{te}=F_{te}^T\hat{\beta }+r_{te}^T{R}_{tr}^{-1}\left(Y^{tr}-F_{tr}\hat{\beta }\right)$$(10)

Where $R_{tr}=\left(R\left(x^{tr}-x{\text{'}}^{tr}\right)\right)$ is the n × n matrix between all training points, $r_{te}^T=\left(R\left(x^{te}-x_1^{tr}\right),\dots ,R\left(x^{te}-x_n^{tr}\right)\right)$ is the vector of correlations between the test point and each training point, F_tr is the n × p known matrix whose i^th row is $\left(f_1\left(x_i^{tr}\right),\dots ,f_p\left(x_i^{tr}\right)\right)$, $1\le i\le n$ and $F_{te}=\left(f_1\left(x_i^{te}\right),\dots ,f_p\left(x_i^{te}\right)\right)$ is the 1 × p vector of known regressors at the test data input.

The regression coefficients and process variance are obtained from:

$$\hat{\beta }={\left(F_{tr}^T{R}_{tr}^{-1}{F}_{tr}\right)}^{-1}{F}_{tr}^T{R}_{tr}^{-1}{Y}^{tr}$$(11)

$$\widehat{{\sigma }_z^2}=\frac{1}{n}{\left(y^{tr}-F_{tr}\hat{\beta }\right)}^T{R}_{tr}^{-1}\left(Y^{tr}-F_{tr}\hat{\beta }\right).$$(12)

The correlation hyperparameters θ are determined by maximizing the log-likelihood of the model:

$$\frac{\partial l\left(\theta \right)}{\partial {\theta }_j}=-\frac{1}{2}\left[Tr\left(R_{tr}^{-1}\frac{\partial R_{tr}}{\partial {\theta }_j}\right)-{\left(Y^{tr}-F_{tr}\hat{\beta }\right)}^T{R}_{tr}^{-1}\frac{\partial R_{tr}}{\partial {\theta }_j}{R}_{tr}^{-1}\left(Y^{tr}-F_{tr}\hat{\beta }\right)\right].$$(13)

Since an analytical solution for θ cannot be obtained directly. An iterative numerical optimization algorithm is employed to minimize the negative log-likelihood. For each candidate value of θ, the corresponding estimates of $\hat{\beta }$ and $\widehat{{\sigma }_z^2}$ are computed using equations (11) and (12), and the likelihood function is updated until convergence.

The Gaussian Process interpolation framework thus provides a statistically grounded and computationally efficient means of approximating aerodynamic responses. Once trained, the GP surrogate model can predict the aerodynamic coefficients or stability derivatives for new configurations almost instantaneously, enabling rapid evaluation of the BWB UAV design space.

4.2 Latin hypercube for sampling

The selection of design points used to construct the surrogate models is known as the sampling process. The accuracy, robustness, and generalization capability of any surrogate model strongly depend on both the number and distribution of the sampled design points within the design space. Since each point corresponds to a computational run of the aerodynamic solver, the sampling strategy must achieve a balance between the representativeness of the data and the overall computational cost. In this context, DoE techniques are commonly used to determine the most effective way to distribute sample points across multiple input variables.

Traditional sampling methods such as full-factorial or fractional-factorial designs systematically evaluate all possible combinations of design variables. While exhaustive, these approaches quickly become infeasible as the number of variables increases, since the required simulations grow exponentially 2^α or 3^α), and the points are typically concentrated at the corners of a hypercube (Fig. 5) [34]. For higher-dimensional problems, LHS [35] provides a more efficient alternative. LHS divides each variable’s range into equally probable intervals and ensures that each interval is sampled once, leading to a balanced and uniform coverage of the input domain. A common guideline is to select a sample size approximately five to ten times the number of variables (N_s = 5 – 10 × α) [35], which provides sufficient data for surrogate construction while keeping computational cost low.

The maximin variant of LHS enhances the overall distribution by maximizing the minimum distance between sample points, thereby improving space filling characteristics and reducing clustering. These properties are particularly advantageous for Gaussian Process (Kriging) models, where uniformly distributed samples lead to more accurate correlation based predictions. By ensuring thorough coverage of the multidimensional design space, LHS enables dependable interpolation and allows the surrogate model to capture both global and local variations in the aerodynamic response.

Fig. 5 Comparison of space-filling properties between Latin Hypercube Sampling and three-level full-factorial sampling in a three-dimensional design space.

4.3 Cross validation

To evaluate the performance of the surrogate models, cross-validation [36] is commonly employed as a robust and reliable statistical method. It assesses the model’s predictive capability on unseen data by dividing the available dataset into two subsets: one used for training the model and the other for validation. This procedure provides an objective measure of the model’s generalization ability and helps prevent overfitting. A key strength of cross-validation lies in its ability to leverage multiple splits of the dataset. There are $\left(\begin{array}{c}n\\ n_p\end{array}\right)$ possible ways to divide the dataset, where n is the total number of data points and n_p denotes the size of the validation set. Several cross-validation strategies exist, such as leave-one-out (LOO) method and k-fold cross-validation.

In the k-fold cross-validation, the dataset is split into k equally sized groups (folds). For each iteration, one group (denoted as T_j) is set aside for validation, while the remaining k - 1 groups are used to train the model. This process is repeated k times, ensuring that each fold serves once as the validation set. The overall predictive performance is then obtained by averaging the errors from all k runs.

The prediction error estimated by cross-validation can be quantified using the Mean Squared Error (MSE), defined as:

$$CV=\frac{1}{n}{\displaystyle \sum }_{j=1}^k{\displaystyle \sum }_{s\in Tj}{\left(Y_s-{\hat{Y}}_s\right)}^2$$(14)

Where s denotes a data point belonging to T_j, Y_s is the simulation value of the target variable, and ${\hat{Y}}_s$ is the model prediction at the same point s.

5 Model selection

In order to select the best model for predicting static stability coefficients a comparison is made to compare the performance of GP (Eq. (7)) with P(x) = β₀, corresponding to ordinary Kriging. The analysis examined the influence of different Latin Hypercube Sampling (LHS) sizes on model performance. Model accuracy was evaluated using a 10-fold cross-validation approach, selected based on prior validation studies [37]. The dataset consisted of seven input variables (Tab. 2), and the corresponding cross-validation results are summarized in Tables 3 and 4.

To quantify the computational benefit of the proposed framework, the wall-clock time of direct XFLR5 evaluations and trained GP predictions was measured on the same workstation. A single XFLR5 analysis required approximately 3 min per configuration, so generation of the 50 points training database required about 150 min. Once trained, the GP model predicted a new configuration in approximately 3 s, which is several orders of magnitude lower than rerunning the aerodynamic solver for each candidate geometry. This confirms that the main computational advantage of the framework lies not only in the use of XFLR5 instead of higher-fidelity CFD, but also in replacing repeated direct aerodynamic evaluations during design screening by near instantaneous surrogate predictions.

Tables 3 and 4 report the 10-fold cross-validation errors obtained for the GP surrogate models with a constant mean function, using the seven-variable parameterization and three LHS sizes, namely 35, 50, and 70 points. Overall, the results show that the 50-point sampling plan yields the best predictive performance for both the longitudinal and lateral-directional stability derivatives. For the longitudinal outputs, the error in C_mα is about 90% higher for the 35-point design and 204% higher for the 70-point design relative to the 50-point case, which confirms that the intermediate sampling level is the most suitable for this response. By contrast, the variation in C_m0 remains very small, with increases of only 0.17% and 9.93% for 35 and 70 points, respectively, indicating a limited sensitivity of this coefficient to the sampling size. A similar behavior is observed for the lateral-directional derivatives where for C_lβ the prediction error is 21.74% higher at 35 points and 13.04% higher at 70 points, whereas for C_nβ the corresponding increases are 17.18% and 29.30%. Taken together, these results indicate that an LHS size of 50 provides the most favorable compromise between prediction accuracy and computational cost among the sampling levels considered.

The relative influence of the design variables on the predicted stability derivatives can then be interpreted through the GP hyperparameters θ, listed in Table 5, which define the correlation structure of the surrogate model through equation (9). In this context, larger values of θ Indicate that the corresponding output is more sensitive to variations in that design variable, and therefore that the variable has a stronger effect on the aerodynamic response [38,39]. Based on the hyperparameter values, Table 6 summarizes the relative influence of each design variable on the stability derivatives.

In summary, the analysis of the Gaussian Process correlation parameters provides valuable insights into the aerodynamic effects of each geometric feature. The root chord and twist difference predominantly govern longitudinal stability behavior, the outward dihedral primarily influences lateral roll stability, and the kink chord exerts a critical effect on yaw dynamics. These findings emphasize the importance of carefully selecting and refining these variables during the design process of a BWB UAV to enhance static stability performance.

6 BWB UAV feasible designs

Based on the cross-validation analysis, the model constructed using 50 LHS points was selected to predict all aerodynamic stability derivatives. As outlined in Section 2, securing static stability for a BWB UAV requires satisfaction of four specific conditions $\left(C_{m_0}=0, C_{m_{\alpha }}<0,C_{l_{\beta }}<0 and{C}_{n_{\beta }}>0\right)$. The flowchart (Fig. 6), illustrates the methodology adopted for design exploration. The search for stable configurations began by pinpointing the zeros of C_m0 within its GP model. We then cross-reference these specific roots within the GP model of $C_{m_{\alpha }}$, $C_{l_{\beta }}$ and $C_{n_{\beta }}$ to ensure compliance with the other stability condition.

The process identified several UAV configurations that satisfied all stability criteria. From these, the three configurations shown in Figure 7 and Table 7 were selected as representative feasible solutions, illustrating the range of stable layouts found within the explored design space. Figure 8 further presents the variation of the three moment coefficients for these configurations, all evaluated about the center of gravity, which remained fixed throughout the analysis, thereby confirming compliance with the prescribed stability requirements. Beyond meeting the theoretical stability conditions, these configurations also provide practical design references for achieving aerodynamic stability in BWB UAV concepts. Moreover, one of the selected configurations had been examined in the authors’ earlier work [30], where aerodynamic predictions obtained with XFLR5 were in good agreement with wind tunnel measurements, while complementary Froude similarity and CAP analyses indicated promising flying qualities [40].

Fig. 6 Post-processing workflow used to identify statically feasible configurations from the surrogate-model predictions.

Fig. 7 Three representative feasible tailless BWB configurations selected from distinct regions of the static-stability-feasible design space.

Fig. 8 Stability-coefficient variations of the three representative feasible configurations: (a) pitching-moment coefficient versus angle of attack, (b) rolling-moment coefficient versus sideslip angle, and (c) yawing-moment coefficient versus sideslip angle.

7 Conclusion

This work developed a unified framework for the early-stage static-stability assessment of tailless blended-wing-body mini-UAVs by combining XFLR5 aerodynamic analysis, Latin Hypercube Sampling, and Gaussian Process Regression. Rather than treating stability as a late verification step, the proposed approach formulates it as a screening problem at the conceptual-design stage, where many geometric alternatives must be evaluated quickly and interpreted physically. From a scientific standpoint, the study shows that surrogate modeling can capture the main relationships between planform geometry and the principal longitudinal and lateral-directional static-stability derivatives with sufficient fidelity for preliminary design exploration. From an engineering standpoint, the framework provides a practical means of identifying feasible tailless BWB configurations before resorting to more expensive CFD or experimental tests.

The results show that the intermediate sampling plan considered in this study offers the best compromise between predictive accuracy and computational effort, and that the resulting surrogate models are sufficiently accurate to support rapid design-space exploration. More importantly, the analysis clarifies the relative roles of the selected geometric variables in governing stability behavior: root chord and root-tip twist difference were found to be the most influential parameters for longitudinal stability, outward dihedral dominated roll stability, and kink chord had a strong effect on directional stability. The identification of feasible configurations satisfying all adopted static-stability criteria confirms that inherent stability can be achieved in tailless BWB mini-UAVs through appropriate geometric shaping alone, even in the absence of conventional tail surfaces.

The present study nevertheless has several limitations. The aerodynamic database is based on XFLR5, which is appropriate for low-speed attached-flow conceptual analysis but does not fully capture viscous effects, separation, or strongly nonlinear flow phenomena. In addition, the investigation is restricted to static stability around a reference operating condition, while other aspects such as dynamic stability, control effectiveness, and structural constraints were not included. The geometric exploration also remains bounded by the adopted seven-variable parameterization, with fixed airfoil selection and center of gravity assumptions.

These limitations define clear directions for future work. A natural extension is to couple the present framework with dynamic stability and control analyses, then integrate aerodynamic, structural, and mission performance objectives into a broader multidisciplinary optimization setting. In addition, the aerodynamic trends identified here should be further examined over a wider range of configurations and operating conditions, particularly in regimes where nonlinear effects may become significant. In this sense, the proposed methodology should be regarded as a computationally efficient tool for conceptual design and static-stability screening, providing a useful basis for subsequent higher-fidelity investigations.

Funding

This research received no external funding.

Conflicts of interest

The authors have nothing to disclose.

Data availability statement

The data generated and/or analyzed during the current study are available within the manuscript and can be provided by the corresponding author upon request.

Author contribution statement

Conceptualization, M. Hakim. and S. Choukri.; Methodology, M. Hakim and S. Coukri.; Software, M. Hakim.; Validation, M. Hakim and S. Choukri.; Formal Analysis, M. Hakim. and S. Choukri.; Investigation, M. Hakim. and S. Choukri.; Resources, S. Choukri.; Data Curation, M. Hakim. and S. Choukri.; Writing – Original Draft Preparation, M. Hakim. and S. Choukri.; Writing – Review & Editing, M. Hakim. and S. Choukri.; Visualization, M. Hakim. and S. Choukri.; Supervision, S. Choukri. and M.E. Ait Ali.; Project Administration, S. Choukri. and M.E. Ait Ali.

References

All Tables

All Figures

	Fig. 1 Body axis system and motion variables.
In the text

	Fig. 2 Aerodynamic forces and pitching moment about center of gravity.
In the text

	Fig. 3 Aerodynamic forces and pitching moment about center of gravity.
In the text

	Fig. 4 Baseline tailless BWB geometry and definition of the seven design variables used in the parametric study.
In the text

	Fig. 5 Comparison of space-filling properties between Latin Hypercube Sampling and three-level full-factorial sampling in a three-dimensional design space.
In the text

	Fig. 6 Post-processing workflow used to identify statically feasible configurations from the surrogate-model predictions.
In the text

	Fig. 7 Three representative feasible tailless BWB configurations selected from distinct regions of the static-stability-feasible design space.
In the text

	Fig. 8 Stability-coefficient variations of the three representative feasible configurations: (a) pitching-moment coefficient versus angle of attack, (b) rolling-moment coefficient versus sideslip angle, and (c) yawing-moment coefficient versus sideslip angle.
In the text