Coupled Aerostructural Optimization of a Composite Low Reynolds Wing Using Surrogate Modeling Techniques

1. Introduction

The wing plays a fundamental role in determining the overall performance of an aircraft, as aerodynamic efficiency and structural weight are strongly coupled through the wing geometry and internal layout. Key performance metrics, such as the lift-to-drag ratio ($L/D$) and take-off weight, are directly influenced by aerodynamic characteristics and structural efficiency. Consequently, careful wing design and optimization are essential to achieving an effective and well-balanced aircraft configuration. During the preliminary design phase, the aerodynamic efficiency of the wing largely dictates the expected performance of the complete aircraft, making early-stage wing optimization a critical step in guiding subsequent design decisions.

Aircraft designers therefore aim to identify optimal configurations from both aerodynamic and structural perspectives. Aerodynamic performance is influenced by parameters such as aspect ratio (AR), taper ratio ($\lambda$), sweep angle ($\mathrm{\Lambda }$), and twist ($ϵ$), while the structural configuration is governed by material selection, internal layout (e.g., ribs and spars), and the thickness of structural components, all of which significantly affect the final weight and structural performance.

Over the past decades, numerous studies have investigated various aspects of wing optimization. From an aerodynamic perspective, Lyu et al. [1] performed aerodynamic shape optimization of a benchmark wing using a gradient-based algorithm coupled with Reynolds-Averaged Navier–Stokes (RANS) equations and the Spalart–Allmaras turbulence model. Similarly, Chen et al. [2] optimized the aerodynamic shape of the Common Research Model (CRM) wing-body-tail configuration under trim constraints. Additional contributions include the work of Ghafoorian et al. [3], who optimized wind turbine blades, and Zheng et al. [4], who applied manifold learning techniques to aerodynamic shape design optimization.

On the structural side, optimization has traditionally focused on minimizing mass while satisfying strength, stiffness, and stability constraints. The introduction of composite materials has significantly expanded the design space by incorporating additional variables such as ply orientation, thickness, and stacking sequence, motivating the development of advanced optimization methodologies [5,6]. Early studies employed direct ply-based optimization approaches, while later work introduced lamination parameters to enable continuous design spaces and improve numerical efficiency [7,8,9]. More recent research has emphasized the importance of geometric nonlinearities in high-aspect-ratio wing structures, demonstrating that linear models tend to overestimate loads and lead to conservative designs [10,11]. Incorporating nonlinear structural behavior has been shown to yield lighter configurations without compromising aerodynamic performance.

Combined, these approaches form aerostructural optimization, which seeks to simultaneously improve aerodynamic and structural performance by accounting for the strong coupling between aerodynamic loads and structural deformation, particularly in flexible wing configurations. Early aerostructural optimization frameworks relied on low-fidelity aerodynamic models, such as lifting-line theory, panel methods, and the doublet lattice method, coupled with simplified structural representations, enabling efficient design space exploration at low computational cost [12,13,14,15,16,17]. With increasing computational capabilities, high-fidelity approaches have emerged, combining Euler or RANS-based CFD solvers with detailed Finite Element Method (FEM) structural models, allowing accurate prediction of aeroelastic effects and load redistribution [18,19,20,21,22]. However, the high computational cost of such approaches limits their applicability during early design stages. Consequently, recent research has focused on multi-fidelity and surrogate-based strategies, which balance computational efficiency and accuracy by combining low- and high-fidelity models [23,24,25].

Surrogate-based optimization (SBO) methods have been widely applied to aerostructural problems. Nikolaou et al. [26] applied SBO techniques to optimize UAV winglet geometry, while Benaouali and Kachel [27] developed a Multidisciplinary Design Optimization (MDO) framework integrating commercial tools, initially focusing on airfoil optimization and subsequently on overall wing performance.

Despite this extensive body of work, the low-Reynolds-number regime remains comparatively underexplored, as most studies focus on conventional subsonic configurations. Furthermore, high-fidelity aerostructural optimization frameworks are computationally expensive and are therefore rarely applied in early-stage design for UAV applications.

To overcome these limitations, the present study proposes a coupled aerostructural optimization framework tailored to low-Reynolds-number UAV wings. Building upon previous work in which a multi-fidelity optimization framework was used to identify an aerodynamically optimal configuration [28], the current study extends this approach by incorporating both aerodynamic and structural considerations, ensuring consistency and comparability of results. The main contributions of this study can be summarized as follows:

• Development of a high-fidelity aerostructural optimization framework integrating CFD and FEM analyses within a surrogate-based optimization (SBO) approach.
• Simultaneous optimization of aerodynamic and structural design variables, including planform parameters and internal wing layout characteristics.
• Application of Kriging surrogate models with an Expected Improvement (EI) strategy to efficiently explore the coupled design space.
• Systematic assessment of surrogate accuracy and optimization robustness through multiple independent runs and validation metrics.
• Identification of the dominant design drivers governing range performance and constraint satisfaction in low-Reynolds-number UAV wings.

The proposed optimization framework employs a high-fidelity surrogate-based optimization (SBO) approach, coupled with Computational Fluid Dynamics (CFD) and Finite Element Method (FEM) analyses. The objective is to maximize aircraft range, subject to constraints on the cruise lift coefficient ($C_{L_{cruise}}$), structural strength via the maximum Failure Index ($F{I}_{max}$), and stiffness via the first global buckling eigenvalue (${\lambda }_1$). The design space includes twelve variables, comprising four aerodynamic parameters—aspect ratio, taper ratio, sweep angle, and tip twist—and eight structural parameters related to internal layout and component thicknesses. Surrogate model accuracy is evaluated using Root Mean Square Error (RMSE) and Leave-One-Out cross-validation metrics, while robustness is assessed through five different Design of Experiments (DoE) strategies.

The remainder of this paper is organized as follows. Section 2 presents the methodological framework, including the CFD and FEM models, the surrogate modeling approach, and the optimization setup. Section 3 reports and analyzes the results of the aerostructural optimization, including surrogate model accuracy and optimal configurations. Section 4 discusses the key findings and concludes the study.

2. Materials and Methods

2.1. Baseline Aircraft Specifications and Requirements

The design process begins with the definition of the mission flight characteristics and the operational requirements of the aircraft. Within the framework of this study, an Unmanned Aerial Vehicle (UAV) has been selected as the case study. According to the NATO classification system for UAVs, the aircraft in this study falls into the Class I mini UAV category. Table 1 presents the key UAV parameters that serve as the basis for the aerodynamic optimization and design framework. The developed UAV is presented in Figure 1.

2.2. High-Fidelity CFD Aerodynamics

A numerical model of a wing, incorporating the average values of key geometric characteristics, was developed as the foundation for subsequent Computational Fluid Dynamics (CFD) analyses and optimization studies of various wing configurations derived from the surrogate model. The computational domain dimensions and boundary conditions were carefully designed to reflect the UAV’s operational environment and altitude. Simulations were conducted using ANSYS Fluent [29], solving the Reynolds-Averaged Navier-Stokes (RANS) equations coupled with the Spalart-Allmaras turbulence model [30]. The Spalart–Allmaras model was selected due to its robustness and proven performance in predicting attached and mildly separated flows at low Reynolds numbers, making it suitable for UAV aerodynamic analysis. The computational domain was defined as a rectangular region measuring 6.0 × 4.0 × 10.0 m. The RANS equations were discretized using the Finite Volume Method (FVM) under incompressible, steady-state flow assumptions with an appropriately refined mesh. To ensure accurate boundary layer resolution, a first cell wall distance of $Y^{+}\approx 1$ was achieved, with the initial layer height set to $y=7.4\times {10}^{-5}$ m. A mesh independence study was also conducted to verify that further refinement did not impact the results, ensuring an optimal balance between accuracy and computational efficiency. The computational domain, including dimensions, boundary conditions, and mesh details, is illustrated in Figure 2 and Figure 3.

Pressure and temperature values were set according to the UAV’s maximum operating altitude, with a predefined inlet velocity. The outlet boundary conditions were defined with a zero pressure gradient, while the turbulence intensity was set to 1%. A symmetry boundary condition was applied along the longitudinal plane, and the wing surfaces were modeled as free-slip walls. The convergence of the lift and drag coefficients, $C_L$ and $C_D$, respectively, is illustrated in Figure 4, while Figure 5 illustrates the $Y^{+}$ distribution for each model of the mesh independence study.

2.3. FEM Model

The FE model of the wing is generated based on the external surfaces of the UAV. The model is generated in NASTRAN, and for the skins, spars, and ribs, 4-noded quadrilateral shell elements (CQUAD4) are used, while for the spar caps and stringers beam elements (CBEAM), accounting also for the relevant offset values. As a datum design point, a cross-ply, two-layered laminate, consisting of a $0^0$ and a ${90}^0$ ply, has been considered for all the relevant wing parts. The Hexcel IM7/8552 composite material system is selected for the relevant parts of the wing, with the respective B-Basis material properties, strength values and cured ply thickness listed in Table 2. For the beam elements and to simplify the analysis, rectangular cross-sections are considered. The respective thickness, nevertheless, is calculated based on the aforementioned laminate, since they are also considered to be manufactured of the provided composite material. However, only isotropic materials are allowed for the definition of the CBEAM elements. As a result, and since the baseline lay-up is symmetric, equivalent laminate axial and shear moduli, $E_{eq}$ and $G_{eq}$ respectively, can be calculated based on the following Equations 1 and [31]:

$$\begin{array}{c}\hfill E_{eq}={\displaystyle \frac{1}{t}}(A_{11}-{\displaystyle \frac{A_{12}^2}{A_{22}}})\end{array}$$

(1)

$$\begin{array}{c}\hfill G_{eq}={\displaystyle \frac{A_{66}}{t}}\end{array}$$

(2)

where t the thickness, $A_{11}$, $A_{12}$, $A_2$ and $A_{66}$ the corresponding terms of the extensional stiffness matrix of a laminate. Regarding the boundary conditions, the wing is assumed to be clamped at its root section, thus fixing all relative nodal Degrees of Freedom (D.o.F). The resulting FEM mesh of the wing model is presented in Figure 6.

2.4. Surrogate Modeling

The key steps in a typical SBO process, as outlined in Alexandrov et al. [33], include:

• Sampling the design space and evaluating the objective function along with any constraints.
• Constructing the surrogate model based on the sampled data.
• Searching the design space and refining the surrogate model using update (infill) criteria.
• Enhancing the model by incorporating newly added points and repeating the process.

The sampling stage is a crucial step in an SBO algorithm, as the surrogate model’s accuracy depends on the selection of initial design points. To ensure the model represents the design space effectively, the most influential points must be chosen to maximize the information available for surrogate construction. Given the often high-dimensional nature of design problems, exhaustive grid searches become computationally prohibitive. Instead, more efficient techniques, such as Latin Hypercube Sampling (LHS) [34], are commonly used. LHS is a robust statistical method that generates parameter samples from a multidimensional distribution while maintaining a well-distributed design space representation. The method involves an optimization problem aimed at maximizing the distance between sample points while ensuring each coordinate follows a predefined probability distribution. Once the sampling is completed, the next step is to construct the surrogate model, typically represented as a general function:

$$\widehat{f}(\mathbf{x},w)$$

(3)

where w denotes model parameters, and $\mathbf{x}$ represents the design variables. A key criterion for selecting a surrogate model is its ability to accurately capture the desired function’s characteristics while maintaining flexibility. Overly rigid models risk instability and overfitting. One widely used surrogate model in engineering applications is Kriging [35,36], which expresses function approximations as a linear combination of basis functions (kernels) that depend on the Euclidean distance between design points. For noise-free data, the Kriging approximation is given by:

$$\widehat{f}\left(\mathbf{x}\right)=\sum _{i=1}^N{w}_i\psi \left(\parallel \mathbf{x}-{\mathbf{x}}^{\left(i\right)}\parallel \right)$$

(4)

where:

• $N_c$ is the number of basis functions,
• ${\mathbf{x}}_c^{\left(n\right)}$ represents the center of the n-th basis function,
• $\psi \left(\right||\mathbf{x}-{\mathbf{x}}_c^{\left(n\right)}|\left|\right)$ is the kernel function, evaluated based on the distance between the prediction point $\mathbf{x}$ and the corresponding center.

The kernel function is typically defined as:

$$\psi \left(\parallel \mathbf{x}-{\mathbf{x}}^{\left(i\right)}\parallel \right)=exp\left(-\sum _{k=1}^d{\theta }_k{\left|x_k-x_k^{\left(i\right)}\right|}^{p_k}\right)$$

(5)

where ${\theta }_n$ and $p_n$ are model parameters.

The Kriging model is constructed using the following steps:

• Formulating the correlation matrix based on training data points:

  $${\left[\mathsf{\Psi }\right]}_{ij}=exp\left(-\sum _{k=1}^d{\theta }_k{\left|x_k^{\left(i\right)}-x_k^{\left(j\right)}\right|}^{p_k}\right)$$

  (6)
• Maximizing the Maximum Likelihood Estimator (MLE):

  $$ln\left(MLE\right)=-{\displaystyle \frac{n}{2}}ln\left({\widehat{\sigma }}^2\right)-{\displaystyle \frac{1}{2}}ln\left(\left[\mathsf{\Psi }\right]\right)$$

  (7)

  where $\widehat{\sigma }$ is the MLE estimate of the standard deviation.
• Predicting values at new design points:

  $$\widehat{y}\left(\mathbf{x}\right)=\widehat{\mu }+{\mathbf{r}}^T\left(\mathbf{x}\right){\mathsf{\Psi }}^{-1}(\mathbf{y}-\mathbf{1}\widehat{\mu })$$

  (8)

where $\mathbf{r}\left(\mathbf{x}\right)$ is the correlation vector between the prediction point and the sampled data points. For Gaussian-based processes, the Mean Squared Error (MSE) estimation is given by:

$$\widehat{s}{\left(\mathbf{x}\right)}^2={\widehat{\sigma }}^2\left(1-{\psi }^T{\left[\mathsf{\Psi }\right]}^{-1}\psi +{\displaystyle \frac{1-{\mathbf{1}}^T{\left[\mathsf{\Psi }\right]}^{-1}\psi }{{\mathbf{1}}^T{\left[\mathsf{\Psi }\right]}^{-1}\mathbf{1}}}\right)$$

(9)

A commonly used approach for improvement is the Expected Improvement (EI) function:

$$E\left[I\left(\mathbf{x}\right)\right]=(y_{min}-\widehat{y}\left(\mathbf{x}\right))\mathsf{\Phi }\left({\displaystyle \frac{y_{min}-\widehat{y}\left(\mathbf{x}\right)}{\widehat{s}\left(\mathbf{x}\right)}}\right)+\widehat{s}\left(\mathbf{x}\right)\varphi \left({\displaystyle \frac{y_{min}-\widehat{y}\left(\mathbf{x}\right)}{\widehat{s}\left(\mathbf{x}\right)}}\right)$$

(10)

where $\mathsf{\Phi }$ and $\varphi$ denote the cumulative distribution and probability density functions, respectively. As an additional step toward more realistic SBO frameworks, constraints should be incorporated into the surrogate model. The approach to constraint handling depends on the computational cost of evaluating the constraint function. Constraints can either be evaluated directly or modeled using surrogate techniques similar to those applied to the objective function, effectively creating a surrogate model for each constraint. When constraint evaluations are computationally inexpensive, conventional constraint optimization methods, in conjunction with the objective function surrogate model, guide the SBO framework toward both promising and feasible regions of the design space. However, if surrogate models are also employed for constraints, the expected improvement function from Equation (10) is modified into the constrained expected improvement function by introducing the probability of feasibility:

$$P\left[F\left(\mathbf{x}\right)\right]=\mathsf{\Phi }\left({\displaystyle \frac{0-\widehat{g}\left(\mathbf{x}\right)}{{\widehat{s}}_g\left(\mathbf{x}\right)}}\right)$$

(11)

where F represents the feasibility measure of a constraint g, and $\widehat{s_g}$ denotes the variance of the constraint’s Kriging model. The probability of achieving an improvement over the current minimum function value while satisfying feasibility conditions is then determined by multiplying Equations (10) and (11):

$$E\left[I\right(\mathbf{x})\cap F(\mathbf{x}\left)\right]=E\left[I\right(\mathbf{x}\left)\right]·P\left[F\right(\mathbf{x}\left)\right]$$

(12)

To determine the next point for model refinement, a sub-optimization problem is solved:

$${\mathbf{x}}_{\mathrm{infill}}=argmax(E\left[I\left(\mathbf{x}\right)\right]·P\left[F\left(\mathbf{x}\right)\right])$$

(13)

This function depends on the Kriging model parameters, ${\theta }_n$ and $p_n$; while optimizing both the correlation length parameters, ${\theta }_n$ and the exponents $p_n$ can improve prediction accuracy in many applications, and doing so increases the complexity of the optimization process. To simplify the model calibration and reduce computational costs, we follow the approach commonly adopted in the surrogate modeling literature [35,36,37] and fix the value of $p_n=2$, corresponding to a Gaussian correlation function. This allows us to focus on tuning ${\theta }_n$ alone, using a global optimization method, which, in our case, is a genetic algorithm, as implemented in Python. The particular optimization was executed for 1000 iterations. As noted by Forrester et al. [37], searching for ${\theta }_n$ on a logarithmic scale, typically within the bounds ${10}^{-3}$ to ${10}^2$, is effective. It is also recommended to scale the input design space to [0, 1] to ensure consistent interpretability of the parameter values across different problems.

The resulting point, ${\mathbf{x}}_{\mathrm{infill}}$, is then added in the current dataset, which is then re-trained. This process is typically repeated for a predefined number of iterations. Within the present study, a samples-to-infill points ratio of 1:2 was selected as recommended in [36]. This ensures that the surrogate model is updated efficiently, balancing exploration and exploitation for improved optimization performance.

2.5. High-Fidelity Aerostructural SBO Framework

The SBO framework employs high-fidelity analysis tools and is executed using five distinct DoE strategies. Following a previous study by the authors, multiple wing configurations are generated using the selected optimized airfoil [28]. The objective of the SBO framework is to identify the optimal aerostructural wing geometry that maximizes the range R, (Equation 14) of the specific aircraft, subject to constraints on the cruise lift coefficient, static strength and stiffness (global buckling).

$$R={\displaystyle \frac{3.6}{g}}·{\displaystyle \frac{L}{D}}·{\displaystyle \frac{E_{sb}{\eta }_{b2s}{\eta }_p{m}_b}{m}}$$

(14)

where:

$$\begin{array}{cc}\hfill m_b& =\mathrm{mass}\mathrm{of}\mathrm{batteries}\left[\mathrm{kg}\right]\hfill \\ \hfill m& =\mathrm{aircraft}\mathrm{total}\mathrm{mass}\left[\mathrm{kg}\right]\hfill \\ \hfill E_{sb}& =\mathrm{battery}\mathrm{specific}\mathrm{energy}[\mathrm{Wh}/\mathrm{kg}]\hfill \\ \hfill {\eta }_p& =\mathrm{propeller}\mathrm{efficiency}\hfill \\ \hfill {\eta }_{b2s}& =\mathrm{total}\mathrm{system}\mathrm{efficiency}\mathrm{from}\mathrm{battery}\mathrm{to}\mathrm{motor}\mathrm{output}\mathrm{shaft}\hfill \end{array}$$

The optimization considers twelve key geometric design variables, comprising four aerodynamic parameters—aspect ratio, taper ratio, sweep angle, and tip twist—and eight structural parameters, including skin thickness, rib thickness, front spar thickness, rear spar thickness, rib spacing, stringer spacing, front spar location, and rear spar location. Each DoE consists of 120 wing models generated within the defined design space. Based on these samples, the SBO framework is trained and subsequently used to generate an additional 60 candidate wing configurations, with the objective of identifying the most suitable design that satisfies the imposed constraints while maximizing range. At the conclusion of each SBO iteration, an optimized wing configuration is obtained, and the resulting optimal designs from the five DoEs are finally compared to assess consistency and robustness of the optimization results. At each SBO iteration the framework starts from the generation of the geometry of the wing, the subsequent computational domain and CFD mesh. The CFD analyses are then conducted, and output is collected in terms of the lift constraint as well as the pressure field in the skins of the wing. Moving to the FEM modeling framework, interpolation schemes are used to map the pressure field from the CFD to the FEM mesh. Given the geometric parameters of the wing, Patran Command Language (PCL) scripts are used to generate a new wing geometry and mesh. Subsequently, analyses files (linear static and global buckling) are extracted and executed, followed by output collection in terms of strength and stiffness constraints. The overall flowchart of the framework is presented in Figure 7.

A summary of all variables, including their lower and upper bounds, is given in Table 3. A horizontal line is used to separate aerodynamic from structural variables. The overall optimization problem setup is summarized in Table 4.

The SBO framework, as illustrated in Figure 8, is common to the previous framework and consists of two main stages: the sampling stage and the model updating stage. The process begins with the definition of the sampling size, followed by generating samples using the LHS method. Subsequently, the geometry and mesh of the geometry are generated and further analyzed via each computational tool. The objective and constraint functions are then obtained for each sample. Once the training stage is completed, the main SBO framework is initiated. The hyperparameters of the Kriging model representing the objective and constraint functions are determined. Next, the constrained expected improvement function (Equation (12)) is minimized using a sub-optimization routine, yielding a new point in the design space. The computational analysis is then performed for this new point, and the surrogate model is updated accordingly. This iterative process continues for a predefined number of infill points.

3. Results

3.1. Surrogate Model Accuracy

The predictive accuracy of the surrogate models was assessed using the Root Mean Square Error (RMSE) and the coefficient of determination ($R^2$). These metrics were evaluated using the infill points generated during the optimization process, which serve as an independent validation dataset. They quantify the average prediction error and the overall agreement between the predicted (${\widehat{y}}_i$) and high-fidelity ($y_i$) responses, and are defined as

$$\mathrm{RMSE}=\sqrt{{\displaystyle \frac{1}{N}}\sum _{i=1}^N{\left(y_i-{\widehat{y}}_i\right)}^2}$$

(15)

$$R^2=1-{\displaystyle \frac{{\sum }_{i=1}^N{\left(y_i-{\widehat{y}}_i\right)}^2}{{\sum }_{i=1}^N{\left(y_i-\overline{y}\right)}^2}}$$

(16)

where N denotes the number of validation points and $\overline{y}$ is the mean of the corresponding high-fidelity responses. A low RMSE and an $R^2$ value approaching unity indicate high predictive accuracy of the surrogate model.

The results for five independent SBO runs are summarized in Table 5.

Overall, the surrogate models achieved a mean RMSE of 5.62 with a standard deviation of $1.11$. For four out of five runs (Runs 2–5), the coefficient of determination remains consistently high ($R^2>0.90$), indicating a strong correlation between the surrogate predictions and the corresponding high-fidelity responses.

In contrast, Run 1 exhibits a significantly lower $R^2$ value, despite having an RMSE of comparable magnitude to the other runs. This behavior is attributed to the relatively narrow range of objective function values in the corresponding infill dataset, which reduces the variance of the reference data and increases the sensitivity of the $R^2$ metric. As a result, even moderate absolute prediction errors lead to a disproportionately low (or negative) $R^2$ value. This behavior is commonly observed in surrogate-based optimization when validation points are concentrated near optimal regions. Therefore, the RMSE provides a more reliable indicator of predictive performance in this case. Considering the consistently low RMSE values across all runs, the surrogate models are shown to provide an accurate approximation of the high-fidelity response, ensuring reliable performance within the surrogate-based optimization framework.

To further evaluate the generalization capability of the surrogate model, a Leave-One-Out (LOO) cross-validation analysis was performed. In this approach, each sample in the DoE is temporarily excluded from the training set, and the model is reconstructed using the remaining $N-1$ samples. The excluded point is then predicted, and the corresponding prediction error is recorded. This process is repeated for all N samples, providing a robust estimate of model accuracy while reducing the risk of overfitting. The LOO Root Mean Square Error (${\mathrm{RMSE}}_{\mathrm{LOO}}$) is computed as

$${\mathrm{RMSE}}_{\mathrm{LOO}}=\sqrt{{\displaystyle \frac{1}{N}}\sum _{i=1}^N{\left(y_i-{\widehat{y}}_{i,\mathrm{LOO}}\right)}^2}$$

(17)

where ${\widehat{y}}_{i,\mathrm{LOO}}$ denotes the prediction obtained by excluding the i-th sample from the training set. In addition, the corresponding coefficient of determination, $R_{\mathrm{LOO}}^2$, was also evaluated in order to quantify the overall agreement between the cross-validated predictions and the high-fidelity responses.

The computed values for the five independent SBO runs are summarized in Table 6.

The LOO errors are highly consistent across all runs, with a mean ${\mathrm{RMSE}}_{\mathrm{LOO}}$ of $8.31\times {10}^0$ and a mean $R_{\mathrm{LOO}}^2$ of $0.818\pm 0.006$. The relatively small standard deviation of both metrics indicates stable surrogate behavior across the independent SBO runs. Although some isolated samples exhibit larger local deviations, as reflected by the maximum absolute errors, the overall cross-validation performance confirms that the Kriging surrogate provides a reliable approximation of the high-fidelity response throughout the sampled design space.

3.2. Best Configurations of All Five SBO Iterations

In this subsection, the best feasible configurations from each SBO iteration are presented and compared to each other. Figure 9 illustrates the convergence history of the cumulative best feasible range (R) for the five SBO iterations. Among the runs, 5^th SBO iteration achieves the highest final range ($\simeq 202.8$), followed by 1^st iteration ($\simeq 196.8$), whereas 2^nd, 3^rd and 4^th iterations converge to similar intermediate values ($\simeq 189-191$).

Table 7 presents the optimal configurations obtained at each SBO iteration, along with their corresponding aerodynamic and structural design parameters. All five configurations exhibit high lift-to-drag ratios ($L/D>25.5$), while the weight remains nearly constant across all cases ($W\simeq 0.83$), indicating that performance improvements are primarily driven by aerodynamic efficiency. Regarding the constraints, $g_1$ and $g_4$ are negative and close to zero, confirming that they are active and govern the feasible boundary of the optimization. In contrast, $g_2$ and $g_3$ remain significantly negative, indicating that they are inactive and do not influence the optimal solutions. From an aerodynamic design perspective, all configurations are characterized by relatively high aspect ratios, with the fifth design achieving the largest value ($AR=14.86$), approaching the upper bound of the design space. The taper ratio and sweep angle remain moderate ( $\lambda <0.44$ and $\mathrm{\Lambda }<{10}^{\circ }$), while tip twist varies between $-0.{40}^{\circ }$ and $2.{11}^{\circ }$, reflecting adjustments in load distribution and aerodynamic performance. In terms of structural design, the parameters are generally consistent across all configurations, with the most notable variations observed in rib spacing (RS) and stringer spacing (SS), suggesting that these variables play a key role in accommodating the aerodynamic changes while maintaining structural feasibility.

Figure 10 illustrates the planform geometry of the five best wing configurations, including their key dimensions: semi-span ($b/2$), root chord ($C_{root}$), and tip chord ($C_{tip}$), but also the ribs, spars (grey lines) and stringers (grey dashed lines). Figure 11 presents the corresponding aerodynamic performance obtained from CFD analyses, including the variations of lift coefficient ($C_L$), drag coefficient ($C_D$), moment coefficient ($C_m$), and lift-to-drag ratio ($L/D$) with angle of attack (AoA). Among the configurations, the 5th wing demonstrates the best overall aerodynamic performance. It achieves the highest lift-to-drag ratio and maintains superior behavior across the entire AoA range, indicating improved aerodynamic efficiency and stability characteristics compared to the other designs.

Therefore, the 5^th wing configuration is selected as the optimal design, as it demonstrates superior performance among the five candidates. It achieves the maximum range of $\simeq 203km$ while also exhibiting the best overall aerodynamic characteristics. In addition, it satisfies all the constraints imposed in the SBO optimization process.

3.3. Selected Optimized Configuration

In this subsection, the results of the 5^th SBO iteration are presented, from which the best optimized wing configuration was selected. In addition, the results of the aerodynamic and structural analyses of the optimized wing are presented. Figure 12 presents the correlation matrix for the final dataset after incorporating all samples and infill points, providing a consolidated view of the relationships among the performance variables and constraints. As observed in all previous iterations (Figure A1, Figure A5, Figure A9 and Figure A13), L/D and R maintain a nearly perfect positive correlation (0.99), confirming that aerodynamic efficiency consistently dominates range performance throughout the optimization process. The correlation between weight and the performance variables remains weakly negative, with W showing correlations of- 0.05 with L/D and -0.18 with R, indicating that, within the explored design space, variations in weight have a relatively limited influence on range compared to aerodynamic efficiency. The constraint correlations remain generally weak with respect to the primary performance variables. Constraint $g_1$ shows almost no correlation with L/D(0.01) and R(0.04), while maintaining a modest negative correlation with weight (-0.25), suggesting limited coupling with the main design variables. Constraint $g_2$ exhibits slightly stronger positive correlations with L/D(0.30) and R(0.33), indicating a mild dependence on aerodynamic performance, while $g_3$ continues to show weak correlations with all variables. Among the constraints, $g_4$ remains the most strongly influenced by the design variables, displaying a moderate negative correlation with weight (-0.65) and moderate positive correlations with R(0.29) and $g_2$(0.64). Overall, the final correlation structure confirms that aerodynamic efficiency is the primary driver of range, while most constraints remain weakly coupled to the performance variables, with $g_4$ showing the most noticeable dependency within the constraint set.

Figure 13 presents the correlation between wing design variables and range R for the fifth and final iteration, reflecting the fully converged design space. Consistent with all previous iterations (Figure A2, Figure A6, Figure A10 and Figure A14), the aspect ratio remains the dominant driver of range, exhibiting a very strong positive correlation ($\simeq 0.94$). The tip twist also maintains a moderate positive correlation ($\simeq 0.43$), reinforcing its role as the most influential secondary aerodynamic variable. In contrast, the taper ratio continues to show a moderate negative correlation ($\simeq -0.30$), while the sweep angle retains only a weak negative influence, confirming its limited impact on range. The influence of structural variables becomes more clearly defined but remains relatively secondary compared to aerodynamic parameters. Stringer spacing shows a small positive correlation, while rib spacing and front spar location exhibit only weak positive effects. On the other hand, several structural variables display consistent negative correlations with range, most notably rear spar location ($\simeq -0.22$), skin thickness ($\simeq -0.27$), rear spar thickness ($\simeq -0.10$), and rib thickness ($\simeq -0.14$). These trends suggest that increases in structural thickness and certain placement parameters tend to reduce range, likely due to associated weight penalties. Overall, the final iteration confirms a stable and well-defined relationship structure: aerodynamic variables—particularly aspect ratio and tip twist—dominate range performance, while structural variables exert smaller, mostly negative influences. This indicates that, at convergence, the optimization has clearly identified the primary performance drivers and reduced uncertainty in the role of secondary design variables.

Figure 14 shows the correlation of aerodynamic and structural design variables with L/D(a) and weight W(b) for the 5^th SBO iteration, reflecting the fully converged relationships in the design space. For L/D(Figure 14-a), the aspect ratio remains the dominant parameter, exhibiting a consistently strong positive correlation ($\simeq 0.9$5). The tip twist retains a moderate positive correlation ($\simeq 0.4$), confirming its role as the most influential secondary aerodynamic variable. The taper ratio continues to show a moderate negative correlation, while the sweep angle has only a weak negative effect. Compared to earlier iterations, the influence of structural variables on L/D is minimal. Most structural parameters cluster close to zero, with only small negative correlations observed for rear spar location, skin thickness, and spar thicknesses, indicating a weak detrimental effect on aerodynamic efficiency. For weight W(Figure 14-b), structural variables clearly dominate. Skin thickness maintains a very strong positive correlation ($\simeq 0.85$), confirming it as the primary driver of weight. Other structural variables, such as front spar location, rib spacing, and stringer spacing, show small negative correlations, while spar thicknesses and rib thickness exhibit weak positive or near-zero effects. Aerodynamic variables have negligible influence on weight, with correlations remaining close to zero or weakly negative. Overall, the final iteration demonstrates a well-converged and decoupled relationship structure, where L/D is governed almost entirely by aerodynamic variables (primarily aspect ratio and tip twist), and weight is dominated by structural thickness—especially skin thickness—with minimal cross-coupling between aerodynamic and structural design variables.

Figure 15 presents the correlation of aerodynamic and structural design variables with the constraint functions $g_1$, $g_2$, $g_3$ and $g_4$ for the 5^th SBO iteration, reflecting the fully converged constraint–design relationships. As in previous iterations (Figure A4, Figure A8, Figure A12 and Figure A16), all constraints are defined such that $g_i<0$ corresponds to feasible designs; thus, negative correlations indicate improved constraint satisfaction, while positive correlations indicate a tendency toward violation.

For $g_1$(Figure 15-a), tip twist remains the dominant variable, exhibiting a strong positive correlation, confirming that increasing twist drives the design toward violating this constraint. Aspect ratio shows a weak negative correlation, indicating a slight improvement in feasibility with increasing aspect ratio. Structural variables exhibit minimal influence, with only very small correlations, suggesting that $g_1$ is primarily governed by aerodynamic variables at convergence.

For $g_2$(Figure 15-b), a clearer aerodynamic influence emerges. Aspect ratio, taper ratio, and sweep angle show moderate positive correlations, indicating that increasing these parameters tends to reduce feasibility. In contrast, structural variables show mixed but generally weak effects, with front spar location exhibiting a noticeable negative correlation, suggesting a potential role in improving constraint satisfaction.

For $g_3$(Figure 15-c), the correlations remain relatively weak overall. Aerodynamic variables show small positive correlations, while structural variables display mixed behavior. Rear spar location and front spar thickness show positive correlations, whereas skin thickness and rib thickness show negative correlations, indicating that increasing structural thickness helps satisfy this constraint.

For $g_4$(Figure 15-d), the strongest and most consistent relationships are observed. Aerodynamic variables—particularly aspect ratio, taper ratio, and sweep angle—show moderate positive correlations, indicating that improvements in aerodynamic performance tend to push the design toward constraint violation. In contrast, skin thickness exhibits a strong negative correlation, with additional negative contributions from spar thicknesses and rib thickness, confirming that structural sizing remains the primary mechanism for restoring feasibility.

Overall, the final iteration demonstrates a well-defined and decoupled constraint structure: $g_1$ and $g_2$ are primarily influenced by aerodynamic variables, $g_3$ shows weak and mixed sensitivity, and $g_4$ remains the most critical constraint, governed by a balance between aerodynamic drivers (causing violation) and structural thickness variables (ensuring feasibility).

Table 8 and Table 9 show the feasible designs of the 5^th SBO iteration, for both samples and infills, respectively. Constraints $g_1$ and $g_4$ remain close to zero for the feasible configurations, indicating that they primarily govern the optimal solutions. In contrast, $g_2$ and $g_3$ exhibit significantly negative values across all configurations, suggesting that they are inactive and do not restrict the optimization. The most feasible designs are characterized by by a high lift-to-drag ratio values and a relatively narrow range of weight, confirming that improvements in aerodynamic efficiency drive the range optimization. The infill results further demonstrate convergence towards a well-defined region of the design space. Overall, this indicates that the optimization successfully identifies feasible, high-performance configurations, where the trade-off between maximizing range and satisfying the most critical constraints - $g_1$ and $g_4$ - defines the optimal design boundaries.

3.3.1. Aerostructural Analysis of the Best Configuration

This subsection presents the aerodynamic analysis results obtained at an angle of attack of ${12}^{\circ }$, which were used as the basis for the FEM analyses. Figure 16 illustrates the contours of the $Y^{+}$ distribution over the upper and lower surfaces of the selected wing configuration, where the values range from 0 to 3.24.

Figure 17 presents the pressure coefficient and temperature contours on the upper wing surface obtained from the aerodynamic analysis at ${12}^{\circ }$ AoA. A minor flow separation is observed near the trailing edge, along with a small vortex forming at the wing tip.

The previously discussed flow detachment is more clearly demonstrated in Figure 18, which presents the airflow pathlines for the examined case. In Figure 18(b) and (d), a small vortex forms near the trailing edge, characterized by flow rotation about the y-axis, resulting in localized flow separation. Additionally, Figure 18-(c) provides a clearer view of the tip vortex.

The structural response of the optimal wing configuration is also examined, providing important insight in the structural behavior of the optimized design and confirm the effectiveness of the proposed aerostructural optimization framework. In particular, Figure 19 illustrates the displacement field of the wing under aerodynamic loading. The deformation pattern is dominated by bending, with maximum deflection occurring at the wing tip, which is consistent with typical cantilever wing behavior. The magnitude of deformation remains within acceptable limits, ensuring that aerodynamic performance is not significantly degraded due to excessive aeroelastic effects. This result further demonstrates that the optimized design achieves a suitable compromise between structural flexibility and stiffness.

The distribution of the maximum failure index across the wing structure is also presented in Figure 20. The highest values are observed near the wing root region, where bending moments are largest due to aerodynamic loading. Despite these localized peaks, the maximum failure index remains below the allowable limit, confirming that the optimized structure satisfies the strength constraints. A similar trend is observed for the maximum beam stress distribution within the structural members (Figure 21). Elevated stress levels are again concentrated near the wing root and along primary load-carrying components such as spars, reflecting the load transfer mechanism within the wing. The stress distribution is smooth and does not exhibit abrupt concentrations, suggesting that the structural layout and sizing variables have been appropriately tuned during the optimization process. The absence of excessive stress peaks further confirms that the design achieves an efficient balance between weight minimization and structural integrity. Overall, the structural response of the optimal configuration demonstrates that the proposed SBO framework successfully captures the interaction between aerodynamic loading and structural behavior. The results confirm that structural variables, particularly skin thickness and internal layout parameters, are actively driven by strength and stability constraints, while aerodynamic variables govern performance. This interplay leads to a structurally efficient and aerodynamically optimized wing design.

4. Discussion and Conclusions

This study presented a coupled aerostructural optimization framework for the preliminary design of a low-Reynolds-number composite UAV wing, integrating high-fidelity Computational Fluid Dynamics (CFD) and Finite Element Method (FEM) analyses within a surrogate-based optimization (SBO) approach. The methodology enabled efficient exploration of a multidisciplinary design space, combining aerodynamic planform variables with structural sizing parameters, while significantly reducing the computational cost associated with repeated high-fidelity simulations.

The results demonstrated that the proposed SBO framework is capable of consistently converging toward feasible high-performance solutions across multiple independent optimization runs. The surrogate models exhibited strong predictive capability, as confirmed by RMSE and Leave-One-Out validation metrics, ensuring reliable approximation of the underlying high-fidelity responses throughout the optimization process.

From a design perspective, the results revealed a clear separation of roles among the design variables. Aerodynamic parameters, particularly aspect ratio and tip twist, were identified as the primary drivers of range performance due to their direct influence on lift-to-drag ratio. In contrast, structural variables—most notably skin thickness—played a critical role in satisfying strength and buckling constraints, thereby defining the feasible design space. The optimal solutions emerged from the interaction between performance-driven aerodynamic variables and constraint-driven structural sizing, highlighting the importance of properly accounting for aerostructural coupling even in preliminary design stages.

The optimal wing configuration achieved a maximum range of approximately 203 km while satisfying all aerodynamic and structural constraints, demonstrating the effectiveness of the proposed framework for early-stage UAV design. Furthermore, the consistency of the results across different Design of Experiments (DoE) strategies indicates robustness with respect to initial sampling, reinforcing the reliability of the SBO approach.

Despite these promising results, several limitations remain. The present study relies on steady RANS simulations and linear structural analysis, which may not fully capture unsteady aerodynamic effects or geometric nonlinearities in highly flexible wings. Additionally, uncertainties related to material properties, manufacturing tolerances, and operational conditions were not considered in the optimization process.

Future work will focus on extending the proposed framework to address these limitations. In particular, the incorporation of geometrically nonlinear structural models and unsteady aerodynamic simulations would enable more accurate prediction of aeroelastic behavior. The integration of uncertainty quantification and reliability-based design optimization (RBDO) techniques represents another important direction, allowing for more robust and realistic design solutions. Furthermore, the application of multi-fidelity strategies could further improve computational efficiency by combining low- and high-fidelity models within the SBO framework. Finally, the extension of the methodology to full aircraft configurations, including fuselage and tail interactions, would provide a more comprehensive assessment of overall aircraft performance.

Overall, the proposed aerostructural SBO framework provides a robust and efficient tool for the preliminary design of low-Reynolds-number UAV wings, bridging the gap between high-fidelity analysis and computationally tractable optimization.

Appendix A.1. 1^st Iteration

Figure A1. 1^st SBO iter. - Correlation between objective and constraints

Figure A1. 1^st SBO iter. - Correlation between objective and constraints

Figure A2. 1^st SBO iter. - Correlation of aerodynamic and structural design parameters with Range-Objective

Figure A2. 1^st SBO iter. - Correlation of aerodynamic and structural design parameters with Range-Objective

Figure A3. 1^st SBO iter. - Correlation of aerodynamic and structural design parameters with $L/D$ (a), and $Weight$ (b).

Figure A3. 1^st SBO iter. - Correlation of aerodynamic and structural design parameters with $L/D$ (a), and $Weight$ (b).

Figure A4. 1^st SBO iter. - Correlation of aerodynamic and structural design parameters with $g_1$ (a), $g_2$ (b), $g_3$ (c), and $g_4$ (d).

Figure A4. 1^st SBO iter. - Correlation of aerodynamic and structural design parameters with $g_1$ (a), $g_2$ (b), $g_3$ (c), and $g_4$ (d).

Appendix A.2. 2^nd DoE

Figure A5. 2^nd SBO iter. - Correlation between objective and constraints.

Figure A5. 2^nd SBO iter. - Correlation between objective and constraints.

Figure A6. 2^nd SBO iter. - Correlation of aerodynamic and structural design parameters with Range-Objective.

Figure A6. 2^nd SBO iter. - Correlation of aerodynamic and structural design parameters with Range-Objective.

Figure A7. 2^nd SBO iter. - Correlation of aerodynamic and structural design parameters with $L/D$ (a), and $Weight$ (b).

Figure A7. 2^nd SBO iter. - Correlation of aerodynamic and structural design parameters with $L/D$ (a), and $Weight$ (b).

Figure A8. 2^nd SBO iter. - Correlation of aerodynamic and structural design parameters with $g_1$ (a), $g_2$ (b), $g_3$ (c), and $g_4$ (d).

Figure A8. 2^nd SBO iter. - Correlation of aerodynamic and structural design parameters with $g_1$ (a), $g_2$ (b), $g_3$ (c), and $g_4$ (d).

Appendix A.3. 3^rd DoE

Figure A9. 3^rd SBO iter. - Correlation between objective and constraints.

Figure A9. 3^rd SBO iter. - Correlation between objective and constraints.

Figure A10. 3^rd SBO iter. - Correlation of aerodynamic and structural design parameters with Range-Objective.

Figure A10. 3^rd SBO iter. - Correlation of aerodynamic and structural design parameters with Range-Objective.

Figure A11. 3^rd SBO iter. - Correlation of aerodynamic and structural design parameters with $L/D$ (a), and $Weight$ (b).

Figure A11. 3^rd SBO iter. - Correlation of aerodynamic and structural design parameters with $L/D$ (a), and $Weight$ (b).

Figure A12. 3^rd SBO iter. - Correlation of aerodynamic and structural design parameters with $g_1$ (a), $g_2$ (b), $g_3$ (c), and $g_4$ (d).

Figure A12. 3^rd SBO iter. - Correlation of aerodynamic and structural design parameters with $g_1$ (a), $g_2$ (b), $g_3$ (c), and $g_4$ (d).

Appendix A.4. 4^th DoE

Figure A13. 4^th SBO iter. - Correlation between objective and constraints.

Figure A13. 4^th SBO iter. - Correlation between objective and constraints.

Figure A14. 4^th SBO iter. - Correlation of aerodynamic and structural design parameters with Range-Objective.

Figure A14. 4^th SBO iter. - Correlation of aerodynamic and structural design parameters with Range-Objective.

Figure A15. 4^th SBO iter. - Correlation of aerodynamic and structural design parameters with $L/D$ (a), and $Weight$ (b).

Figure A15. 4^th SBO iter. - Correlation of aerodynamic and structural design parameters with $L/D$ (a), and $Weight$ (b).

Figure A16. 4^th SBO iter. - Correlation of aerodynamic and structural design parameters with $g_1$ (a), $g_2$ (b), $g_3$ (c), and $g_4$ (d).

Figure A16. 4^th SBO iter. - Correlation of aerodynamic and structural design parameters with $g_1$ (a), $g_2$ (b), $g_3$ (c), and $g_4$ (d).