Integration aspects of the collaborative aero-structural design of an unmanned aerial vehicle

Overall aircraft design is a complex multidisciplinary process, which requires knowledge from many different fields such as structures, aerodynamics, systems and propulsion. For unconventional configurations lacking an empirical knowledge base, higher fidelity physics-based methods are required to reliably estimate the feasibility of a given new design concept. Analysis tools and results are provided by highly specialized groups of experts, possibly from different organizations. In the AGILE (aircraft 3rd generation MDO for innovative collaboration of heterogeneous teams of experts) project, new approaches to setting up cross-organizational collaborative aircraft design optimization workflows have been investigated, including the employment of common parametric aircraft configuration schema as a central common data schema and the provision of disciplinary analysis competences as callable services. Following this paradigm, the present paper details a distributed workflow to perform an aero-structural design optimization of an unmanned aerial vehicle (UAV) design. Taking advantage of disciplinary capabilities provided by several partners based in various locations across Europe, an integrated design workflow including a distributed and tightly coupled aero-structural analysis loop has been assembled using the process integration and design optimization system remote component environment developed at the German Aerospace Center. To enable the necessary load and displacement transfer between non-matching disciplinary meshes, a versatile and lightweight algorithm using radial basis functions has furthermore been implemented. The functionality of the workflow is demonstrated by performing the optimization on the baseline configuration of the UAV.


Introduction
Interdisciplinary collaboration is becoming increasingly important in overall aircraft design. To mitigate entrepreneurial risks, it is important to incorporate as much knowledge as possible during early design stages. This requires close cooperation between disciplinary experts and aircraft designers. To this aim, the tools and knowledge of disciplinary experts from diverse backgrounds are integrated in automatic design workflows, which can even extend beyond company boundaries with increasing regularity. So far, the large amount of time required to set up these workflows has often had a negative impact on their utility, leaving too little time to run large-scale analyses in time-constrained projects. Therefore, the AGILE Horizon 2020 European project [1] aims to employ new methods, to significantly lower the time required to set up and reconfigure multidisciplinary design workflows [2]. Taking advantage of the expertise of many partners beyond Europe, the goal is to perform multidisciplinary design optimization (MDO) on six different configurations.
While most of the designs considered in AGILE are evaluated on a preliminary level, the unmanned aerial vehicle (UAV) is used as a test case, in how far the methods developed in AGILE can be applied in a higher fidelity context. Therefore, an optimization workflow for the UAV has been implemented, which takes into account static aeroelasticity, coupling Euler computational fluid dynamics (CFD) analyses and global finite element method (FEM) analyses. A special challenge here is each analysis being carried out by a different partner in a tightly coupled iterative process, which requires stable and well-defined interfaces, ensuring correct and structured data transfer even after many iterations.

The OPTIMALE configuration
The baseline for the collaborative multidisciplinary design task at hand is the OPTIMALE configuration (cf. MALE: medium altitude long endurance) developed during the German AeroStruct research project [3]. A representation of the design in the common parametric aircraft configuration schema (CPACS) [4,5] has been provided to the AGILE consortium by Airbus Defence and Space. A visualization of the model using the TiGL viewer [6] is shown in Fig. 1.
As can be seen, the OPTIMALE is a conventional low wing configuration with a T-tail and two rear-mounted turbofan engines. Information on engine performance is provided by the Baranov Central Institute of Aviation Motor Development. The system power offtake is taken into account based on an onboard system architecture designed using the process described by [7].
Top-level aircraft requirements missing from the Aero-Struct baseline have been supplemented by Airbus Defence and Space in collaboration with other partners from the AGILE consortium. A selection of key requirements for the aeroelastic analysis is given in Table 1.
Two reference missions have been defined for the OPTI-MALE. One is a transfer mission where the aircraft is relocated between bases. It consists of a takeoff, climb, cruise, descent and landing segment. The other reference mission is a surveillance mission, which is illustrated in Fig. 2. Differently from the transfer mission, a loitering segment is introduced in addition to the cruise segment, where the aircraft remains in the same location at low altitudes to gather surveillance data. Both cruise and loitering altitudes are given in Table 1.

Problem definition
One of the goals of the AGILE project is to set up a distributed aeroelastic shape optimization process for the OPTI-MALE. The optimization problem to be solved can be stated as follows:  The objective function of the optimization is the endurance for the surveillance mission t surv , which must be maximized. Since the mission segments for takeoff, landing and cruise are fixed, this implies a maximization of the potential duration of the loitering segment. During the computation of t surv , the nested optimization problem described in Sect. 3.3.2 is solved to size the airframe structure, adding thickness and strength constraints to the overall optimization problem. Several wing planform parameters are given as design variables, which are illustrated in Fig. 3. As can be seen, the span b, and chord length c i and twist i for three span-wise positions i per wing are modified.

Downselection of flight load cases
Since the deformations of the aircraft and the resulting structural loads are unique to each flight condition, a separate multidisciplinary analysis (MDA) loop must be run for each case. Given the considerable time required to run the CFD analysis especially, this calls for a downselection of the considered load cases. For the analysis of the OPTIMALE, it was decided to focus on what are presumably the two most critical maneuver load cases within the envelope. Consequently, a 3 g and a − 1.5 g maneuver at dive speed v d at sea level were chosen. The corresponding flight conditions including Mach number Ma and lift coefficient C L used in the aerodynamic calculation are listed in Table 2. In addition to the sizing cases, the conditions for loitering at 1.05 g are also given. In the following, the individual building blocks of the workflow will be discussed in detail. Due to the particular architecture used, which features nested loops and parallel calls, the work setting up the present process has been performed manually using the process integration and design optimization system remote component environment (RCE) [9,10], instead of the MDO automation tools in AGILE, which do not support ad hoc MDO architectures yet.

Disciplinary model generation and analysis
Since the workflow execution is distributed among multiple partners, special emphasis must be placed on clearly designing the interfaces, where communication between partners is required. Most importantly, a commonly accepted unambiguous data model of the aircraft design to be investigated is necessary. Here, the CPACS format is a very good choice, since it not only provides a fairly detailed geometrical description of the aircraft, but also supplementary data on masses, mission, engine performance, etc. As CPACS is also an increasingly popular format in aircraft design, automatic model generators are available to produce consistent aerodynamic and structural models from a single source.
The structural model is generated by Airbus Defence and Space using Descartes [11], a versatile interactive tool suite for working with CPACS files developed in-house, which  provides the capability to build global finite element models based on the structural description from a given CPACS file. The models are written as input files for the Lagrange software [12], an FEM solver and optimizer also developed in-house by Airbus Defence and Space. The input format of Lagrange is very similar to the Nastran format in many respects, which means that existing tools for Nastran can be adapted with manageable effort. The tool is not only used to perform the linear-static analyses in the fluid-structure MDA, but also to provide optimization capabilities, which are taken advantage of during the subsequent sizing optimization loop. A key advantage of Lagrange is the availability of analytical derivatives with respect to the design variables, mitigating the need for expensive numeric differentiation. Several optimization algorithms are available; however, all analyses were performed using a sequential quadratic programming algorithm with distributed and non-monotone line search [13].
The aerodynamic analysis components are provided by CFS Engineering and Airinnova. The generation of the aerodynamic surface mesh is accomplished using the sumo software [14]. Since sumo does not accept CPACS as input, a compatible data set must be generated from the baseline using the CPACS2SUMO interface [15]. The volume mesh is grown from the surface mesh using TetGen [16]. The CFD analysis itself is performed using SU2 [17], an open-source CFD solver developed at Stanford University. To reduce the run time of the analysis and the complexity of the mesh, Euler analysis was chosen over RANS. Since the lift coefficient C L for a given flight condition is prescribed according to Table 2, it is necessary to perform the analysis in fixed C L mode, where the angle of attack is varied to receive a given value for C L . SU2 also provides methods for volume mesh deformation, which can propagate displacements on the surface mesh to the volume mesh using the inverse volume stiffness method. As the structure is deformed by the aerodynamic loads, the deformations must be propagated back to the volume mesh for the next CFD run during each iteration.
Furthermore, SU2 only computes a surface pressure distribution. Post-processing is necessary to obtain the forces at the grid points, which can be mapped to the structural model. This can be done by computing the normal force for each surface cell using where A,cell is the area-scaled cell normal and p cell is the pressure on the cell surface. The force at a given point can be computed simply by averaging the cell normal forces of all adjacent cells. Figure 5a, b shows the FEM and CFD meshes, respectively, which have been generated from the initial CPACS geometry.

Interface mapping methods
One of the key problems when solving partitioned multiphysics problems is the transfer of interface data between non-matching meshes. In the case of aeroelasticity, the displacements on the boundary of the structural domain must be mapped to the fluid domain and, in return, the forces on the boundary of the fluid domain must be mapped to the Over the past years, meshless surface reconstruction techniques have become the commonly adopted solution to perform the mapping in fluid-structure interaction problems, due to their dependence only on the mesh point locations, but not the mesh connectivity. Generally speaking, the discrete mapping equation computes the displacements fluid on the fluid domain boundary based on the displacements at the structural domain boundary struct using where is the interpolation or mapping matrix. Due to the principle of virtual work, the forces can be mapped as well conserving the energy by using the transpose of the mapping matrix: Thus, the problem of mapping forces and displacement is reduced to the problem of computing the mapping matrix .
Two common solutions to this have been investigated in detail: the interpolation method based on radial basis functions (RBF) described by [18] or [19] and the moving least squares (MLS) method as described by [20].
The RBF method has been widely adopted and solves a global interpolation problem taking into account all model points at once. The mapping matrix is constructed as follows: with and The matrix contains polynomial parameters, e. g., denotes the evaluation of the radial basis function for the Euclidian distance between any two discrete points in the domain boundaries Γ 1 and Γ 2 . The major drawback of the RBF method is that, for larger meshes, the matrices Φ fl,st and Φ st,st can quickly outgrow the memory available in commonly available computer hardware. If compactly supported functions are used, much can be achieved by exploiting that Φ fl,st is a sparse matrix for compactly supported classes of basis functions; however since st,st is a dense matrix, the number of structural nodes remains a limiting factor.
In contrast, the MLS method performs a local leastsquares approximation of a series of subsets of points. For each point in the aerodynamic set, a subset ̄ of the structural point set is picked based on a support radius or a n nearest neighbor criterion. The surface is reconstructed for the subset using a least-squares approximation: where remains the polynomial matrix, is the polynomial parameter vector for a single point, and diag is the RBF evaluation of the diagonalized distance vector between and ̄ i. e., Concatenating ( ) for all in the evaluation set yields the coupling matrix .

Fig. 5 Disciplinary models from CPACS
As it happens, the matrix st,st need not be computed or inverted, which is an advantage over the RBF algorithm. Furthermore, the large global problem is broken down into many smaller and independent local problems, which can potentially be solved in parallel. However, it must be noted that the MLS algorithm can only provide a good approximation of the original data, whereas the RBF method is interpolating in that it reproduces the same data at the same point.
Both algorithms have been implemented in Python by German Aerospace Center (DLR) as part of this project. Due to the higher accuracy, it is, however, preferable to use the RBF method, if overfitting due to large differences in displacement within a small radius is not an issue and the models are sufficiently small. This is the case for the present models. A combination of a Wendland C0 radial basis function and second degree polynomial parameters proved to give the best results and was therefore selected as default coupling method.

Collaborative multidisciplinary analysis
The collaborative MDO workflow is centered around the MDA loop, where the aeroelastic equilibrium between structural deformation and aerodynamic forces is computed in an iterative process employing a Gauss-Seidel fixed-point iteration scheme on the structural deformation.
An initial set of deformations {0} FEM of the structural mesh is mapped onto the aerodynamic surface mesh points using Eq. (3). Conveniently, it is sufficient to compute the mapping matrix once at the beginning of the aeroelastic sizing process, since only the displacements are adapted, while the underlying meshes remain unchanged. The mapped results in comma-separated value format (CSV) are passed to CFS Engineering and Airinnova, who are responsible for performing the aerodynamic analysis, where they are taken as inputs for the volume mesh deformation. To alleviate configuration effects between wing and empennage, only the wing deformations are taken into account so far. The pressure distribution on the surface cells, resulting from the subsequent Euler analysis, is then post-processed into force vectors acting on the individual mesh points as outlined in Sect. 3.1.
The force vector distribution on the aerodynamic surface points is returned to the parent workflow and mapped back onto the structural model using Eq. 4. An updated set of FORCE input cards [21] for the structural solver is written and a structural analysis is performed. The residual is given by the 2-norm of the change in deformation of all FEM nodes. The loop is converged if the following condition is fulfilled: where the absolute tolerance set to = 5e−3. If the condition is not fulfilled, the loop enters another iteration using the current displacements as input. A comparatively high value is chosen as tolerance to keep the analysis time reasonably low. In addition, the number of iterations for the MDA is limited to 8. If the loop does not converge, a structural sizing is performed anyway in an attempt to generate a more fitting stiffness distribution, which may result in more stable deformations.
As mentioned previously, the MDA loop must be executed separately for each flight load case, since each case requires a different CFD analysis to be run and different forces and displacements to be mapped. However, all loops can be run in parallel, to keep the time penalty low. Upon convergence the final loads for all flight conditions are passed to the structural sizing optimization loop.

Structural sizing optimization
Once all subcase loops have converged or reached the maximum number of iterations, the final force distribution is collected and forwarded to the structural sizing optimization, which is performed using Lagrange. The following optimization problem is solved: where the sizing variables cells are the skin thicknesses of the individual buckling fields on the wing, which are sized against strength criteria. The buckling fields not only include the upper and lower cover of the wing, but also the leading and trailing edge skin and the ribs. Since these components consist of isotropic and anisotropic materials, criteria based on both stress and strain are applied. The objective is to minimize the structural mass m struct .

Convergence of the full sizing loop
Following the structural sizing optimization, the deformation of the wing under the given loads, as well as the structural mass of the aircraft are likely to have changed. In this case, it is necessary to reenter the aero-structural sizing loop by returning to the aero-structural MDA. Once the change in mass is sufficiently small, the sizing loop is converged and the aero-structural sizing has completed. A change in mass Δm ≤ 10 kg is considered sufficient for convergence.

Data exchange
As outlined in Sect. 3.1, the CPACS format has been adopted as the common data structure for exchanging aircraft data. However, for the data transfer within the aeroelastic loop most of the capabilities of CPACS cannot be taken advantage of, since exchanges between simulation grids have not been defined. Due to the comparatively large quantities of homogeneous tabular data being exchanged, the XML format furthermore adds unnecessary overhead. Instead, it is preferable to use a dedicated tabular format such as CSV, which can easily be parsed by common numerical libraries. The structure of the tables must be agreed upon by the partners in advance. As an example, the format of a table for displacement exchange is given in Table 3. The format for force transfer is analogous. Python interfaces have been implemented by the disciplinary experts to modify their analyses based on incoming CSV data as well as provide tables as outputs.
The exchange of data across company bounds via RCE is enabled by the BRICS protocol [22,23]. Developed at the Netherlands Aerospace Center, BRICS provides the capability to share files on a neutral server, which can be accessed by all partners in connection with a request system. The calling workflow has to send a work request to the service provider, who needs to accept the request by launching the service manually. This way, the control over if and when the service should be run resides with the provider and unwanted data extraction is prohibited. Originally developed with CPACS in mind, BRICS can also be used to exchange CSV files if the communicating processes are designed accordingly.

DOE/design optimization loop
The aero-structural sizing loop is nested within the aircraft design loop, where the global design parameters are modified and the performance of the aircraft, i.e., the objective function, is evaluated. Before the aeroelastic sizing is performed, the geometry of the baseline design is morphed according to the design variables. The aerodynamic and structural meshes are deformed accordingly, as well.
Once the structural mass has been determined in the sizing loop, the mass breakdown in the CPACS file for the deformed design is updated and a mission simulation is performed.

Shape deformation
The deformations due to the changes in the shape parameters are applied to the CPACS geometry using Descartes. The tool also provides the capacity to update the structural model by deforming the baseline mesh. The surface deformations of the morphed structural models with respect to the baseline are computed by forming the difference between the two node sets and and mapped to the CFD surface mesh using one of the mapping algorithms introduced in Sect. 3.2. For both, the structural and the aerodynamic mesh, the deformations are applied only to the nodes, keeping the mesh topology unchanged. This procedure is beneficial, not only because the meshes need not be regenerated for each deformation, saving run time, but also because it retains the possibility to pass gradient data throughout the multidisciplinary system. On the other hand, the dependency on existing meshes limits the design space, since excessively large deformations will result in poorly conditioned meshes, which may give incorrect results. However, in the scope of this analysis, it is assumed that the initial design is already reasonably close to an optimal solution, so small variations are sufficient.

Performance evaluation
The performance of the design on the surveillance mission is evaluated at DLR using the Aircraft Performance Program (APP) [24,25]. A Python interface has been implemented to extract the relevant masses, engine data and aerodynamic performance data from CPACS and provide them to APP.

Mapping
As a first step, the results of the interface mapping procedure described in Sect. 3.2 are given. Figure 6a shows the forces at the mesh points computed from the pressure distribution of a CFD run on the undeformed aerodynamic mesh. The forces are mapped to the FEM mesh using Eq. (4), resulting in the distribution shown in Fig. 6b. As can be seen, the magnitude is increased, however, since there are fewer application points, the resultant force remains equal.
For the next step, a static analysis, using the force distribution given in Fig. 6b is performed. The resulting

Multidisciplinary analysis results
With the mapping results validated, an isolated run of the aeroelastic MDA is performed for the baseline. Figure 8 illustrates the convergence of the step norm ‖ ‖ Δ FEM ‖ ‖2 over the course of eight iterations for both the 3 g and the − 1.5 g case. While the − 1.5 g case converges one step before the 3 g case, it is shown that the step norm is decreasing continuously for both. The loop converges approximately linearly, especially toward the end, with a convergence rate of approximately 0.3. The convergence criterion given in Eq. (11) is fulfilled in the final step, and the iteration breaks off before the maximum number is reached.

Structural sizing
Based on the converged loads for the 3 g and − 1.5 g cases computed in the MDA in Sect. 4.2, a structural optimization is performed using Lagrange. Figure 9 shows the initial thickness distribution of the model on the left wing and the sizing result on the right. For CFRP shells, the overall laminate thickness is given. As can be seen, the thickness is increased especially near the wing root, whereas it can be reduced toward the wing tip. This is a plausible result, since the wing bending moments can be expected to be the highest at the root.
The overall structural mass is decreased by 40 kg during the optimization, indicating that the initial distribution of the baseline was already realistic. It must be noted, however, that very restrictive gages have been set for the leading and trailing edge thicknesses, to avoid unrealistically large deformations. This results in the leading and trailing edge being

Aero-structural sizing loop
Running the MDA and the structural sizing sequentially in an iterative loop yields the aero-structural sizing loop. In Fig. 10, the results for both the step norm during the MDA and the structural masses after optimization are given. It can be seen that the iteration converges after only two aeroelastic sizing loop runs, fulfilling the criterion given in Sect. 3.3.3. This once again indicates that the initial mass estimation was already fairly close to the optimum under the given design constraints.
It can also be noted that the MDA converges in both runs, and the iteration is not stopped prematurely due to the maximum function call criterion. However, the initial change in displacement after the first structural optimization run is again quite high. This is due to the fact that the displacements for the first deformations are reset to zero again at the start of the MDA. While the reasons for this are process related, it can be assumed that the second MDA loop could have been accelerated by taking the last deformations of the first loop as starting point.

Mission simulation
The surveillance mission is simulated using APP for the unsized baseline configuration. Using the engine performance map supplemented by CIAM and an aerodynamic performance map generated at DLR by running a panel method analysis with viscous drag correction, the results in Fig. 11 were obtained. The mission simulation shows that the OPTIMALE can stay in the air for almost 24 h, traveling a distance of approximately 8250 km at cruise speed v cruise and loitering altitude h loiter from Table 1. The mission is planned with a 5% fuel reserve, which is reflected in the corresponding plot.
While already providing an estimation of the capabilities of the OPTIMALE, the results of the mission simulation will likely be subject to change, as more accurate methods for generating aerodynamic performance data become available.

Mesh permutations for the design of experiments
To perform the DOE, the baseline geometry must be deformed with the respective disciplinary analysis meshes following accordingly, using the methods outlined in Sect. 3.4.1. Figure 12 shows the planform variations on the FEM mesh and the mapping results on the CFD mesh for two example cases. Since it is possible for two neighboring points (e.g., a fuselage and a wing point) to experience very different deformations, the RBF algorithm is liable to overfitting when mapping said deformations to the CFD surface mesh. In contrast, the MLS algorithm provides a smooth surface approximation, which is why it has been selected for this particular use case.

Fig. 10
Step norm and mass during the iteration

Conclusions
A collaborative multidisciplinary design process for the OPTIMALE configuration, which has been implemented as part of the AGILE project, is presented. The work is shared between Airbus Defence and Space, who provide the model baseline and structural analysis and sizing competences, CFS Engineering and Airinnova for CFD model generation and analysis, and DLR for workflow integration, mesh data mapping and mission simulation. The workflow consists of a design of experiments (DOE) component driving a nested aeroelastic sizing loop. Before entering sizing, the geometry is transformed according to the design variables. Endurance is evaluated in a mission simulation, which is performed based on the new geometry and masses after sizing is completed. The sizing loop itself comprises a multidisciplinary analysis (MDA) loop to compute an aeroelastic equilibrium and a structural sizing optimization, which are run in sequence. Two surface mapping techniques to transfer loads and displacements between the aerodynamic and structural domains are discussed as well.
The results in Sects. 4.1, 4.2 show that the individual components work as expected. Furthermore, a full sizing was performed on the baseline. A strategy for applying the geometry deformations from the DOE is presented in Sect. 4.6. This demonstrates that all the building blocks necessary to perform the DOE and explore the design space for the OPTIMALE are in place. However, while individual cases have already been run, the full DOE is yet to be finished.
Upon completion, the next step is to replace the DOE component with an optimizer to solve the optimization problem stated in Eq. (1). To this end, the endurance computed in the mission simulation must be fed back to the optimizer to compute a new set of design variables. However, it must be noted that an optimization will require a substantially larger amount of workflow evaluations than to the DOE. While gradient-based optimization has the potential to reduce the necessary number of function calls significantly over finite differences, it is difficult to integrate in the present workflow, especially due to the nested optimization using Lagrange. Therefore, it may be beneficial to investigate other MDO architectures, possibly involving surrogate models, to circumvent this issue. The automation tools introduced in AGILE have the potential to drastically simplify the task of regenerating of the workflow, once the challenges posed by the data standards in high-fidelity analysis and ad hoc architectures like the present one are overcome.