A framework for the bi-level optimization of a generic transport aircraft fuselage using aeroelastic loads

The aeroelastic loads and design processes at the German Aerospace Center, Institute of Aeroelasticity in the framework of multi-disciplinary optimization are constantly evolving. New developments have been made in the in-house model generation tool ModGen, which allow us to create detailed fuselage models for preliminary design. As a part of the subsequent developments to integrate the fuselage structure in our aeroelastic design process, a new framework for optimizing the fuselage structure has been developed. The process is based on a bi-level optimization approach which follows a global–local optimization methodology to simplify a large optimization problem. A sub-structuring procedure is used to define stiffened panels as independent structures for local optimization. The panels are sized with stress and buckling constraints with consideration of several aeroelastic load cases. Furthermore, in this paper, we present a physical sub-structure grouping process which enables reduced number of panel optimizations and saves considerable computational effort with little compromise in the solution accuracy.


Introduction
Multi-disciplinary optimization (MDO) is an indispensable and integral part of the overall aircraft design process. The goal in developing such MDO frameworks is to obtain a good preliminary estimation of the structural dimensions usually with the objective of weight reduction and also to gain assurance toward the structural design. The complexity of the process and the fidelity of the models used vary depending on the design level for which the framework is developed. For example, MDO at the preliminary design level is focused on the global aero-structural performance and works sufficiently well with medium fidelity models. At the German Aerospace Center, Institute of Aeroelasticity (DLR-AE), a comprehensive framework, cpacs-MONA [1], has been developed for the aeroelastic pre-design of fullscale aircraft models. Until now, the focus of this framework has been on the optimization of the wing structures based on aeroelastic loads and constraints. The fuselage structure is therefore modeled as an equivalent beam structure and is sized using analytical methods on a conceptual design level. Recent developments at DLR-AE in the in-house model generation tool ModGen [2] have enabled the possibility of modeling a full finite element model (FEM) of the fuselage structure utilizing the design data defined in the CPACS (Common Parametric Aircraft Configuration Schema) format [6].
To integrate the fuselage in our design process, a Python-based automated framework has been developed at DLR-AE to conduct the structural optimization of the fuselage based on aeroelastic loads. The fuselage structure is primarily constructed of stiffened panels consisting of shell and beam elements. Due to the large number of elements used to generate the finite element (FE) model of the fuselage, an optimization model can potentially consist of thousands of parameters. The application of mathematical programming techniques for designing such structures can be impractical due to the requirement of high computational power. Furthermore, large-scale optimization problems can become infeasible when using standard optimization algorithms [9]. One approach for reducing the large dimensionality of an optimization problem is to identify weak couplings in a group of variables and divide it through the use of local variables [3,8]. A preliminary estimate of the material distribution can be subsequently obtained through the use of low-fidelity methods such as fully stressed design (FSD). The cross-sectional dimensions are then resized in a detailed design step involving local design variables and more accurate mathematical programming techniques. This method is generically termed as the bi-level optimization approach. A comprehensive description of an optimization framework based on the above principle is presented by Grihon et al. [3]. Such formulations have also been applied to fuselage barrel section models built with composite materials [4] where the failure criteria, such as buckling, damage tolerance and post-buckling, are utilized in the local optimization problem. Furthermore, optimization results for a composite fuselage section including local elastic tailoring of the laminate properties are presented in [14]. A multiobjective optimization methodology for a fuselage section including window and door cut-outs has been discussed in [15], where the authors consider minimizing the deformation of the cut-outs as an additional design objective combined with mass minimization. A multi-scale optimization strategy for fuselage structures using genetic algorithm has been discussed in [20,21]. In these articles, notably, the boundary conditions for sub-structure optimizations are defined without simplifying hypotheses through the use of enforced displacements which are interpolated across the edge grid points of the local models. This approach enables a more realistic coupling between a global and a local model. A similar process has also been utilized in sizing composite wings of varying configurations including fatigue and manufacturing constraints [22][23][24]. The application of bi-level optimization does not necessarily ensure computational efficiency and in certain cases, the solution process for a realistic fuselage model can require several days for completion [18]. This necessitates development of efficient optimization techniques. In correlation to this, the use of surrogate models as a replacement for traditional solution methods in the framework of MDO is an area of interest. An assessment of the performance benefits of applying well-known surrogate modeling techniques to an MDO process for preliminary design of wings has been conducted by Paiva et al. [10]. The authors conclude that the approach is worth investigating however, the benefits are not clearly predictable. It was found that the methods begin to lose advantage with increasing complexity and dimensionality of the optimization problem. Nevertheless, variations of surrogate modeling techniques have been successfully implemented in broader MDO frameworks with good results [12]. A numerical sizing approach using FSD is presented in [16] where the considerations for the local buckling stability is incorporated through handbook methods. The utilization of FSD allows obtaining quick estimates of a design, however with limited accuracies.
In this article, new developments in the aeroelastic design process at DLR-AE are discussed. The focus of this work is on the development of an efficient automated framework for sizing fuselage structures using bi-level optimization techniques which can be coupled to our in-house aeroelastic loads computation process for the preliminary sizing of a full-scale fuselage model of a generic transport aircraft. To alleviate the computational costs in the local optimization, a physical grouping procedure is proposed for the fuselage panels. The grouping procedure allows selective sizing of a reduced number of panels to ensure faster solutions with reasonable accuracy. The framework is applied to a full-length fuselage and the skin thickness distribution is computed. Several maneuver and quasi-static gust loads computed at different flight conditions and mass configurations are utilized in the analysis, and the thickness distribution obtained through optimization is presented.
The fuselage barrel of the D150 aircraft model first developed in the DLR project "Virtual Aircraft Multi-Disciplinary Analysis and Design Process" (VAMP) [5] is used for the studies. The format for defining the structural elements of the fuselage in the CPACS dataset is discussed in [26]. The loads applied on the full model are obtained through static aeroelastic analyses in the SOL 144 module [25] of MSC Nastran which is setup within cpacs-MONA. In the first level, a global finite element model (GFEM) is sized using the fully stressed design module (FSD) of MSC Nastran [13]. A linear static analysis is conducted to extract internal forces and stresses in the full FE model. A substructuring procedure is used to generate stiffened panels for local optimizations. Subsequently, the internal forces are utilized for panel optimizations to obtain a global thickness distribution.
The subsequent discussion can be outlined in the following manner: (1) in Sect. 2, a description of the primary structural components in the FE model and the reference D150 aircraft configuration is presented, (2) Sect. 3 describes the process flow in the optimization framework which includes the model development, aeroelastic load computation, structural optimization and the proposed sub-structure grouping procedure, (3) the application of the framework to the D150 aircraft and the corresponding results are discussed in the Sect. 4, (4) and finally the concluding remarks are presented.

Structural components in the model
The finite element model consists of an assembly of the skin panels stiffened by the longitudinal stringers and the circumferential orbit frames. The stringers are essential in carrying the tension and compression loads on the fuselage, whereas the orbit frames maintain the shape of the fuselage and helps in load redistribution [11]. The skin panels sustain the majority of the cabin pressure loads. In Fig. 1, the skeletal structure of a fuselage section, generated using ModGen, is depicted. In addition to the main fuselage barrel section, the keel beam, pressure bulk heads, passenger and cargo floors are also modeled (Fig. 2).
The keel beam lies in the center fuselage region where the wing loads are introduced into the fuselage structure. The center wing box (CWB) is modeled for the wing load transmission into the fuselage.
The FE models of the fuselage generated through ModGen are low-fidelity structural models for early preliminary design. One shell element is formulated between two parallel stringers and two orbit frames. The tool allows mesh refinement depending on the design level of interest. A sequence of line elements is defined along the fuselage length and the cross-sectional mid-points. This axis is termed as the loads reference axis (LRA) and the corresponding nodal points are used for the introduction of the aeroelastic loads into the structure. The LRA is connected to the encompassing structure through rigid body interpolation elements which enables the load distribution into the primary structure.

Reference aircraft model
The optimization framework presented here has been assessed with the use of the reference aircraft model D150, shown in Fig. 3. This aircraft was designed as a mid-range, single-aisle configuration, transport aircraft for 150 passengers. In general, the aircraft design has similarities in characteristics to the Airbus A320 model. Some of the relevant dimensions and design data of the D150 fuselage are listed in the Table 1. The cross-sectional dimension of the center fuselage barrel is depicted in Fig. 4. It is noted that the stringers, which run from the nose to the tail end of the fuselage, are circumferentially not equally spaced. Therefore, a range for the stringer pitch is listed. Similarly, the frames have variations in the interspacing across the fuselage length .    Two types of cross sections have been defined for the stiffening structures: 'Z' shapes are used for the stringer geometries, while 'C' shapes are used for the frame geometries. The initial dimensions of the cross sections are defined to be constant for each element across the circumference and the length of the fuselage. Figure 5 depicts the relevant dimensions (in mm) for definition of the stiffening elements.
The FE model of D150, built in MSC Nastran, has 10,759 node points and over 25,000 elements. Quadrilateral shell elements (CQUAD4, Mindlin-Reissner shell) with 24 degrees of freedom (DOF) per element are predominantly used to model the main fuselage barrel, whereas the triangular shell elements (CTRIA3, Mindlin-Reissner shell) with 18 DOF per element are used as substitutes in areas where insufficient space leads to poor aspect ratio and distortion of the quadrilateral elements. The stiffening components of the structure, such as the orbit frames and the stringers, are modeled using the 1-D beam elements with 6 DOF per node. The connection between the shell and the beam elements inherently exists due to node commonality. The material properties used in the FE model are described in Table 2. The allowable strengths have been listed with an applied factor of safety of 1.5.

Optimization framework
The framework utilizes the concept of super-stiffeners [3] to break down a large optimization problem into several substructure optimizations. In this section, a description of the sequential steps in the setup of the framework is given with reference to Fig. 6.
The first two steps involving the model generation and aeroelastic loads computation are a part of the external tools ModGen and cpacs-MONA, respectively. A description of the model components has been provided in Sect. 2. In cpacs-MONA, the loads are generated using the aeroelastic solver SOL 144 in MSC Nastran. The flight conditions are defined at various Mach numbers, flight speeds and altitudes ranging up to 12.5 km based on the specification of the certification document CS25 [17]. Additionally, five different mass configurations with forward and aft center of gravity (CG) positions are considered. The variation in the CG positions is achieved through appropriate distribution of payload and fuel masses. For a comprehensive description of the load case definitions, aeroelastic load computation and design process in cpacs-MONA, the reader is referred to [1,7].

Load case selection (LCS)
Several aeroelastic load cases are obtainable through the combination of different flight conditions and mass configurations. However, the utilization of all generated load cases enforces redundant computations in the sizing process. This occurs since each sub-structural component is sized using the corresponding critical load cases i.e., the cases which generate the maximum magnitudes of forces and moments. To circumvent unnecessary computations, a load case selection (LCS) process is applied. In this process, a comparison is made between the integrated loads at each fuselage section. The load cases generating the maximum magnitudes of forces f i and moments m i at any j th fuselage section are determined to form a reduced subset of loads comprising F (1)  Thereafter, the reduced loads subset is applied in the sizing process.

Fully stressed design
The FSD process is applied to the global finite element model (GFEM) of the fuselage to get a preliminary estimate of the thickness distribution. The FSD methodology is not designed to give an accurate representation of the required material, but rather is seen as a lower bound estimate [13]. It is based on the premise that each member of the structure should be at its maximum allowable stress under at least one of the design load conditions. Therefore, the design variables are updated iteratively such that the stress levels reach the maximum allowable in each element of an FE model. In practice, this is implemented as discussed in [19]: where i is the iteration number, N is the total number of design variables and t i j is the current value of the j th design variable, i j is the current stress level for the element which corresponds to the j th design variable and al j is the allowable stress. The SOL 200 module in MSC Nastran is used for this purpose. The selected parameters for the sizing process are: (1) shell thicknesses t of the elements (8345 in total) are selected as design variables, (2) von Mises stresses on the upper and lower element surfaces (16,690 in total) are selected as the design responses, (3) and the minimum skin thickness in each iteration is defined to be 1 mm which fulfills the requirement of defining design constraints together with the allowable stresses.
This keel beam and the center wing box are presently only sized using the FSD process and are not considered in the optimization model. Obtaining a realistic initial thickness estimate in the preliminary FSD step is also necessary for accurate internal force computation, which directly influences the subsequent steps in the process flow. An example of the sizing results obtained from the FSD process is shown in Fig. 7, where a thickness gradient can be seen on the keel beam. The center wing box is sized with a single design variable for each panel, therefore a uniform thickness is obtained through this process.

Internal force extraction
A linear static analysis of the fuselage model is conducted in SOL 101 module of MSC Nastran using the design properties obtained from the FSD process. The internal forces and moments generated at every node point in the FE model due to the applied external loads are evaluated in this step. These internal loads are required for the local panel optimizations. The analysis is conducted using the inertia relief feature in MSC Nastran since the fuselage structure is globally unconstrained during operational conditions. The inertia relief feature balances out the external forces and moments by applying nodal accelerations, thus maintain a state of perfect equilibrium for an unconstrained structure. Therefore, no assumptions with respect to definition of boundary conditions are required in the global analysis.

Super-stiffener generation
The super-stiffener element is a stiffened shell panel which has a width equal to the stringer pitch in the GFEM and the length equal to the distance between two orbit frames. Figure 8 depicts a super-stiffener element on a larger fuselage panel and the corresponding extracted super-stiffener geometry.
According to Grihon et al. [3], isolating such a superstiffener geometry to analyze the structure locally is a good approximation when a suitable stiffening in the transverse direction is provided. This concept is well applicable to fuselage structures due to the presence of orbit frames. Through application of the super-stiffener extraction procedure on the main fuselage body, 6593 independent stiffened panel geometries can be created for the D150 aircraft model.

Grouping of super-stiffener geometries
It is evident that conducting optimization for all generated panels would require a high computational time. The grouping process has been implemented to reduce the number of super-stiffener panel optimizations that need to be conducted independently to obtain a fully updated optimized solution. To this end, two solution methods were investigated: loadbased grouping and von Mises stress-based grouping.
All super-stiffener panels are first segregated into groups based on similar dimensions. In the load-based method, the internal forces and moments at the node points of the GFEM are compared. The panels having a similar range of forces and moments are assigned one super-stiffener group. Although, this approach is more generically applicable, it is more tedious due to the requirement of local coordinates for every panel and the necessary comparison of all three components of forces and moments for the appropriate grouping.
The alternate approach using the von Mises stress as the grouping criterion was found to be more effective. In this case, normalized stresses in every super-stiffener panels are compared to formulate panel groups. In every group, the panel which is subjected to the maximum stress is optimized and the new properties are applied to all super-stiffener panels of that group.
This procedure is suitable for structures constructed with isotropic materials for which the optimization is generally conducted with von Mises stress as one of the design constraints. In the case of composite materials, where the sizing is often conducted with principal strain as a design constraint, it is more practical to allocate groups based on similar loads or strains. The fundamental idea in any of the  grouping procedures remains the same while the difference only arises from the choice of the grouping criterion. Figure 9 shows an example of the stress-based grouping. The number of panels lying within a certain stress interval is shown with the blue markers. The enclosing curve defines the maximum number of grouped panels in every stress interval. The choice of stress interval is objectively decided by the user depending on the extent of error margin acceptable in the design. A conservative design is obtained when a high stress interval is chosen as the grouping criterion. Further discussion on this is presented in the Sect. 4.

Local optimization
The local optimization is conducted on the super-stiffener panels, as discussed previously in Sect. 3. The panel optimizations are conducted with constraints on the von Mises stresses such that the stress allowable of the material a is not exceeded. Additionally, it is ensured that the critical load factor , obtained through an eigenvalue buckling analysis of the panels, is greater than 1. The internal forces extracted through the static analysis of the GFEM is utilized in this step. These forces are applied on the edge 1 and edge 3 of the panel shown in Fig. 8b. The internal forces are directly transferable to the grid points that are common in the global and local models. Since the local panel model has a finer meshing than the GFEM, the loads are distributed on to the panel edges using the RBE3 interpolation elements of MSC Nastran. The key characteristic of the RBE3 element is that it does not add additional stiffness to the structure, however, can be used to transfer forces and moments from a reference master node to its dependent nodes. The remaining edges 2 and 4 are simply supported to constrain the structure from undergoing rigid body motion in static analyses. It is notable that assumptions in the boundary conditions do not necessarily reflect the global structural behavior accurately [21]. Furthermore, edge effects may also influence the analysis for which some measures have been recommended in [20,21]. These considerations are currently not incorporated into the framework presented here and therefore, may lead to some deviations from global analyses.
With the defined optimization model, the panels are optimized sequentially and the new properties are then allocated to the global FE model. The process is iterated with the updated properties until convergence in the thickness distribution is obtained. The optimization is conducted in the SOL 200 solver [13] of MSC Nastran which utilizes a gradient-based optimizer. Design sensitivities are computed for the objective function and the design constraints which provides a search direction for the optimization algorithm. A relative change in the objective function of < 0.001 between two design cycles is chosen as the convergence criteria for the algorithm.

Application and results
The D150 aircraft is used in the subsequent studies to demonstrate the applicability of the framework. Initially, one aeroelastic load case from cpacs-MONA labeled as LC1, is considered to size the fuselage. The grouping procedure is used in this case and a comparison of the sizing results and analysis time is presented for this load case. Subsequently, the complete sets of loads computed at different flight conditions, labeled as LCM, are utilized in the sizing process. The sizing results are finally compared to LC1.

Sizing with one load case: LC1
The load case is computed with the following analysis parameters:  The internal pressure load is added as the difference between the cabin pressure and the external atmospheric pressure. A reference optimization result is generated without grouping such that all super-stiffener panels are independently optimized. A linear static analysis prior to the optimization process gives an indication of the stress distribution with the initial thicknesses, as shown in Fig. 10. The thickness distribution of the optimized model is expected to have a similar trend.
The obtained thickness distribution from LC1 is shown in Fig. 11. The highest thicknesses are localized near the center fuselage section with a gradual reduction in thicknesses when moving toward the front or rear end of the fuselage. A vast majority of the element thickness lie in the range of 1-2 mm while the maximum thickness in the lateral region of the center fuselage is closer to 3.6 mm. The highest thicknesses obtained are localized in the keel beam.

Influence of stress interval
The stress interval determines the number of panel optimizations that need to be conducted.
As discussed in Sect. 3, the choice of higher stress interval reduces the number of panel optimizations to be done and therefore, reduces the computational time. However, the approximation that the optimization result of an individual panel with the highest stress is valid for all super-stiffeners in that particular group leads to over-estimation in the design. The grouping process is, to an extent, a trade-off between the computational efficiency and solution accuracy.
To understand the effectiveness of the grouping process, a verification study was conducted with different stress intervals and the resultant weight estimates and thickness distributions were compared. The following Table 3 gives an overview on the different weight estimates obtained with various stress intervals. The reference dataset describes the case where all panels are independently optimized and represents the ideal case.
A considerable reduction is seen in the completion time of the optimization cycle when the super-stiffener panel grouping process is applied. However, a corresponding increase in the estimated weights is also observed. This is explainable by observing the thickness distributions obtained with different stress intervals.
Due to the grouping approximations, thicknesses in some panels are over-estimated and therefore, a higher weight estimation is expected. Figure 12 shows the thickness distributions in the optimized models with normalized stress Fig. 11 Thickness distribution of the D150 fuselage obtained from the optimization framework (top, side and bottom view). No grouping procedure is used and a single 2.5 g maneuver load case LC1 is applied intervals: 0.1, 0.01, and 0.001. These plots are generated with the same color spectrum and therefore, give an accurate indication of deviations in thickness distribution. When using a normalized stress interval of 0.1 for grouping, the solution is very quick but the thicknesses in a larger area near the center fuselage are over-estimated, which results in the 11.26% gain in the weight of the optimized model. At the stress interval of 0.01, the thickness distribution has some visible differences from the reference model, and an error of 1.59% in the weight estimation is seen. However, the solution time reduces by 83.6%. At a stress interval of 0.001, the thickness distribution is much closer to the reference model, the error in weight estimate is 1.37% and a reduction of 57.7% in solution time is achieved. There is clear advantage of using the grouping procedure in terms of computational efficiency while producing an acceptable level of accuracy.

Sizing with multiple load cases: LCM
As previously discussed in Sect. 3, cpacs-MONA [1] generates a diverse set of load cases which are utilized in the aeroelastic design process. A fully assembled configuration of the D150 model, including the structural stiffness of wings and the tail plane, is considered in the aeroelastic analysis in SOL 144 in MSC Nastran. The following summarizes the load case parameters that are utilized in the subsequently presented analysis. Varying combinations of the listed parameters result in 693 load cases for the D150 aircraft. In this analysis, the LCS procedure, previously discussed in Sect. 3, produces 24 sizing relevant load cases. The resultant thickness distribution obtained is presented in Fig. 13. The thickness distributions are plotted using the same color spectrum as in Fig. 11. In the middle section of the fuselage, the thicknesses are perceivably higher in comparison to the thickness distribution obtained with a single load case in LC1. This occurs since several super-stiffener panels are now sized by the critical load cases in the corresponding fuselage sections chosen by the LCS process. A constraint monitoring process is utilized to assess the main driver behind the design in different sections of the fuselage. In Fig. 14, panel locations are depicted for which the design is primarily driven by the buckling constraints. In these panels, buckling load factors computed in the final design cycle within the optimization algorithm are found to be marginally higher than the critical buckling load. These locations are predominantly in the central fuselage region and below the center line of the fuselage, especially in the vicinity of the center wing box cut-out. This is in alignment with the expectation that, due to the higher loading in the center fuselage region, the panels are more prone to buckling. Overall, the necessity of considering a comprehensive set of loading conditions in the sizing process is highlighted through an increase of 5.4% in the normalized weight. The loading conditions and design constraints considered here, while necessary, are not invariably sufficient for a detailed design. Further investigations would be carried out with considerations of additional structural constraints, such as fatigue and additional loading conditions, derived from landing and ground maneuvers. As a final remark, the direct comparison of the thickness distribution results to the reference A320 model cannot be published due to confidential data.

Conclusion
A framework for optimization of the fuselage structure of a generic passenger aircraft has been developed. The method utilizes a multi-level optimization approach where an initial thickness distribution is obtained from a lowfidelity model using fully stressed design. A sub-structuring procedure is used to generate sub-components called super-stiffener panels. These panels are locally optimized with stress and buckling constraints using gradient-based optimization in MSC Nastran. A physical grouping procedure for the super-stiffener panels is implemented to reduce the overall computational effort. It is found that the grouping procedure significantly speeds up the process while the errors in the thickness and weight estimates due to grouping are within acceptable margins. While the framework considers some necessary design parameters, further improvements can be introduced. Investigations pertaining to the local optimizations are one of the focal points. In particular, the appropriate transformation from a global to local models is a major point of interest. In the context of the applied load cases, the influence of the landing and ground loads is yet to be investigated. The application of the framework so far has been limited to the D150 aircraft model. A change in the structural model is not expected to have any impact on the described process.