Finite element method on topology optimization applied to laminate composite of fuselage structure

3.1 Problem formulation

In general topology optimization problem, an objective function F is set, which is subject to a volume constraint G 0 ≤ 0 and possibly N other constraints G i ≤ 0 , i = 1 … M . The material distribution is described by the density variable ρ ( x ) , which has a value of either 0 (void) or 1 (solid material) at any point in the design domain Ω. Mathematically, this optimization problem may be written as follows:

subject to G 0 ( ρ ) = ∫ Ω ρ ( x ) d V − V 0 ≤ 0 ,

where the state field u satisfies a linear or non-linear state equation. In this case, the objective function can be seen as the integral over a local function f ( u ( ρ ) , ρ ) , for instance the strain energy density when considering compliance optimization. Elaborated problem formulation of the strain energy density itself is accessible in the study of Hu et al. [17] and its solution approach in virtual work is discussed in the study of Elvas et al. [18].

Eq. (1) can be approached in two possible ways. The first is a Lagrangian approach which is known for its boundary following mesh, while the latter formulation is the Eulerian approach with its fixed mesh. Both approaches have been studied thoroughly in the study of Tong et al. [19] pointing out a concern in the lack of their solutions. In Lagrangian, generating more holes for optimization solutions may jeopardize the object function. In Eulerian, the increase in several constraints may develop similar hole problems.

Another issue is related to the discrete value of the density variable (0 or 1) prompting difficulties in finding solutions, particularly for the Eulerian method. This issue can be tackled with efficient gradient-based optimization algorithms [20]. The problem can be further seen as a continuous topology optimization problem. The algorithm itself is designed to ensure convergence within a reasonable number of iterations. In this new paradigm, the problem can be written as follows:

subject to G 0 ( ρ ) = ∑ i υ j ρ j − V 0 ≤ 0 .

where ρ denotes the design variable vector of length N. Also, the part f ( u ( ρ ) , ρ ) of Eq. (2) explains the relationship between density and material property that can be written as f ( u ( ρ ) , ρ ) = g ( ρ ) f 0 ( u ) , where g ( ρ ) is the density interpolation function and f 0 ( u ) is a function of the field for solid material. The approach described earlier is called the nested approach in which the equilibrium equations are assumed to be satisfied for each optimization step.

An issue of checkerboard problems still may be encountered, which refer to patches of alternating black and white elements as a result of bad finite element modeling but are favored by the unrestricted optimization process. The checkerboards may be avoided by the use of higher-order elements, or they are taken care of by all the restriction methods that ensure mesh independence [21].

On the condition that the structural constraints are met, the objective function of topology optimization minimizes the structural strain energy. Increasing the structure’s stiffness will lower the structural strain energy. Giving each finite element unit density using design variables will enable this concept. Strain energy (U) formulation can be stated as follows:

For instance, when calculating the maximum static structural stiffness with volume limitations, the maximum structural stiffness entails determining the minimal statistical strain energy for a given load, which translates to determining the minimum compliance.

3.2 Topology optimization framework

Figure 2 shows the flow chart of the proposed optimization method. The whole process is divided into three steps. The first stage is the process of evaluating the structure using FEM based on the structural model before the optimization process is carried out.

Figure 2

Work flow on numerical topology optimization.

The second stage is the optimization process which begins with the extraction of the frame structure from the complete global finite element method (GFEM) model. After that, the design response is determined using the parameter values for strain energy and mass weight.

In the optimization process, constraints are also given by providing the optimum model with a minimum limit value of 0.5% of the original mass weight. If the initial model cannot meet the constraints, then the constraints are adjusted. After that, the optimization process can be run by giving the command a limit of 50 design cycles. It should be noted that topology optimization achieves convergence values in the range of 25–40 design cycles. The optimization process stops until the minimum strain value is obtained until it converges. Finally, the structural assessment phase was carried out again by incorporating several optimized frames into a complete airframe structure.

3.3 MALE UAV load generation according to STANAG 4671

One of the critical inputs in the optimization process is the aircraft load analysis. The design, requirements, and objectives (DRO) govern everything related to the product development process, and all certification requirements are specified in the DRO. Subpart C of STANAG 4671 specifies almost all of the requirements for MALE UAV load. The load on the aircraft is divided into flight load and ground load. Flight load includes any load that occurs during airborne operations, while ground load includes any load that occurs during the take-off and landing phases.

All external loads acting on the aircraft must undergo analysis. STANAG 4671 Subpart C specifies the probability of any load that must be taken into account. The aircraft load must be analyzed using the GFEM. The output of the GFEM process is a load at the component level, which is stated in the nodal load. The scheme of aircraft load generation is shown in Figure 3.

Figure 3

Aircraft load generation diagram.

Performing an identification of the potential loads that occur on the aircraft is crucial to obtain the appropriate aircraft load. A weight and balance analysis of the aircraft is an excellent starting point to obtain comprehensive load information by defining the mass distribution throughout the aircraft. Weight and balance analysis is performed by assuming the aircraft load factor to be equal to the maximum load factor. The aerodynamic load consists of lift and drag, both when flying symmetrically and when performing various maneuvers. The aircraft configuration is also a consideration in carrying out load analysis. It is necessary to consider how the load path of the existing load will affect the loading philosophy as stated in the regulations. Attention must be paid to the redundant or single load path fail-safe design philosophy. Similarly, the design of the fuselage, wing shape, sweepback, taper ratio, and tail design play an important role in the load on the aircraft. It is essential to consider the position of the tail on the fuselage, whether it is a conventional tail or T-tail, or whether there is a canard (horizontal tail on the front fuselage). All of these design aspects result in different load variations on the aircraft.

Landing gear provides a critical load on the aircraft. When the take-off phase, the load on the landing gear is relatively not so large and quite steady. But when landing, the impact that occurs on the aircraft due to the collision of the landing gear with the ground becomes quite large and stated as critical load. On this MALE UAV, the landing gear is all located in the fuselage hull. In general, in aircraft design, the load allocation for the nose landing gear is around 10–15%, while for the main landing gear, it is around 85–90%. So, usually, the failure of the fuselage structure occurs at the nose landing gear attachment when the nose landing is down.

The position of the engine and its thrust also impose significant loads on the aircraft. As an aircraft component, the engine has considerable weight, and thus, its position and connections require extra attention. The engine’s thrust generates a compressive load on the fuselage structure from the rear, while the drag generates a compressive load from the front. Therefore, in general, the fuselage structure is under compression.

The high-lift device plays a crucial role in controlling the aircraft. It includes ailerons and/or flaperons in the wings, rudder in the vertical tail, and elevators in the horizontal tail. The movement of the high-lift device does not directly affect the fuselage structure, as it is restrained by the wings and tail. However, the movement of the high-lift device indirectly causes maneuvers that change the load factor on the aircraft, thereby significantly increasing the load experienced by the fuselage.

The mission profile of the MALE UAV needs to be considered to determine the aircraft load. Factors such as flight distance, trajectory, operating altitude, and mission type will determine how the load will occur. The fit between components on the plane is also important as it serves as a load path for internal loads on the plane. The variation of this fitting causes the load on the aircraft to differ. In general, the load source on the aircraft is shown in Figure 4.

Figure 4

Detail source of aircraft loading.

After identifying all potential loads on the aircraft, it is necessary to establish load limits that may occur on the aircraft. As shown in Figure 5, the flight envelope diagram describes the load limits on the aircraft. The left side is limited by the stall speed, which is the minimum speed required for the aircraft to maintain flight. The upper part shows the positive load factor limit that the aircraft can handle due to maneuvers or gust loads. On the right side, the aircraft speed limit is shown, indicating the maximum speed during cruising or diving. The dive speed occurs when the aircraft cruises at maximum speed and then performs a downward maneuver, resulting in an increase in kinetic energy, which is a change in altitude potential energy. Meanwhile, at the bottom, the flight envelope is the maximum limit for a negative load factor. By defining the flight envelope diagram, the designer can determine the limits of the load on the aircraft that should not occur.

Figure 5

Flight envelope diagram of MALE-UAV.

4.1 Benchmark topology without considering orientations

In this section, we describe a numerical prediction method that was validated using several studies found in the open literature. A benchmark example of the MBB beam was taken for the validation of the present topology studies. For this example, a model composite reference from Lee [16] is demonstrated. The detail of the benchmark model is presented in Figure 6. Optimization problems include minimizing adherence to the total structural volume of beam structure under static load and constraint with residual volume as a design variable.

Figure 6

Benhmark model of MBB beam structure.

The mechanical properties of T300/5208 Carbon material used in calculations are listed as follows: E ₁ = 132 GPa, E ₂ = E ₃ = 10.8 GPa, G ₁₂ = G ₁₃ = 5.65 GPa, G ₂₃ = 3.38 GPa, v ₁₂ = v ₁₃ = 0.24, v ₂₃ = 0.59. The MBB beam has been performed with the 2D element. In the load scenario, the forces are applied downward as concentrated load in the mid-bottom center as shown in Figure 7. Both ends at the edge are set to fixed boundary conditions. The MBB beam has been modeled using 2D elements. Mesh size is set to 2 mm achieving the number of elements of 3,750. In this numerical modeling, the angle of plies is fixed to 0°.

Figure 7

Load and boundary conditions on MBB beam benchmark model structure.

The converged topology optimization took only 27 design cycles as presented in Figure 8. A significant reduction can be achieved as shown in Figure 9. The qualitative validation of the MBB Beam subjected to lateral forces shows a good agreement with the benchmark model. The results show the compliance minimization problem with the fixed angles in the optimization with a 40% mass constraint.

Figure 8

Optimization result of the benchmark model.

Figure 9

A significant reduction can be achieved as shown in Figure 9 from (I) to (IV).

4.2 Mesh convergency analysis

Before performing the main analysis, the convergence test was performed by changing the mesh sizes of the modeling to select the proper mesh size. In this analysis, the boundary conditions are set using fixed conditions in frame nos 8 and 9 in the wing lug area where the metal parts are mounted. At this stage, the weight of the instrument is concentrated evenly on each area where the instrument is placed. The magnitude of the applied load is equal to each instrument’s weight such as battery, FCU, PCU, Payload, etc. The convergence test results are shown in Figure 10. As shown in Figure 11, when the element number was 178,571, the maximum deflection converged to 19 mm, which was almost half of the coarse mesh size 70 mm × 70 mm. From this investigation, the selected element size for the main analysis was 50 mm × 50 mm.

Figure 10

Mesh convergence test.

Figure 11

Detail of mesh size: (a) mesh size 10 mm with 178,571 elements, (b) mesh size 30 mm with 44,887 elements, (c) mesh size 50 mm with 16,875 elements, (d) mesh size 60 mm with 12,946 elements, and (e) mesh size 70 mm with 9,662 elements.

5.1 Model and material

The FE model of the fuselage structure has a total length of 8 m with 13 frames and one fire wall at the end of the fuselage construction. Details of the structure are shown in Figure 12.

Figure 12

Schematic illustration of MALE UAV structure configuration.

The material properties of the fuselage structure are summarized in Table 2. The composite laminate is designed using material from HexPly W3G282-F593 carbon fiber prepreg with the composition of several plies with stacking sequence described in Table 3.

The procedure in material modeling is carried out with the following stages. Initially, several layers in the composite section are modeled according to the direction of the layers and the number of layers in the section as shown in Table 3. All composite layers are arranged to form a balanced layer and a quasi-isotropic material. Subsequently, an equivalent-single material of the final composite component is taken to be applied as the homogeneous material in the final analysis. This method is proven to produce robust analysis and reduce computational costs.

5.2 Loads on fuselage

This article presents the loads on the MALE UAV fuselage as follows. The loads on the fuselage are assumed to be applied to each frame along the fuselage. The MALE UAV fuselage consists of 14 frames. Figure 13 shows a typical load for each frame, which is classified into several categories, such as aerodynamics, inertia, attachment force, and fitting force.

Figure 13

Typical load for frame structure evaluation.

These loads are generated by several factors such as mass and distribution, aerodynamics, aircraft geometry, landing gear, engine and thrust, stability and control, mission profile, and fitting. To obtain representative results in topology optimization, the critical flight load and critical ground load have been taken into consideration (Table 4).

Aerodynamic load, inertial load, and attachment force are then represented by a single-point load, which is located in the floor beam of each frame. This single-point load is subjected to several loads and moments and distributes the load to certain points in the frame. Meanwhile, the fitting force is calculated based on the landing gear load data.

6.1 Topology optimization results

In this study, topology optimization required up to 30 iteration steps to meet the convergence criterion. Figure 14 shows an example of how the iteration progress optimized the strain energy and volume fraction of frame No. 1, along with an illustration of the evolving structure. In this example, after 28 iteration steps, the strain energy had been reduced to 17 J from an initial value of 128 J, while the volume fraction had decreased to 67% from its initial value. The decrease in volume fraction was associated with weight reduction.

Figure 14

Strain energy and volume fraction vs number of iterations.

Figure 14 gives an interesting outlook that the convergence of strain energy and volume fraction occurred in different steps. In this particular instance, the volume fraction and strain energy decreased right away as the optimization was initiated. However, the volume fraction has not changed ten iteration steps earlier than strain energy reaching convergence. Here, it can be concluded that the optimization process focused on accomplishing constraint criterion first before settling for the solution of the objective function. Small changes in the topology have shown a significant impact on the strain energy. In structural analysis, strain energy has a physical meaning as the energy stored in a body due to deformation. In this case, the end form is a frame with some holes. It can be seen that these holes transformed from hooked to more curved. The optimization direction of curve-shaped holes is expected as they indicate less stress concentration; hence, less deformation propagation is likely to develop.

As previously mentioned, the decrease in volume fraction associates with the weight reduction; yet, the relationship is not quite linear. In frame no. 1, a 33% reduction in volume fraction gives a 42% weight reduction. Similar occurrences are present in other frames as shown in Table 6. This finding is expected as the density of the frames is not completely uniform; thus, the decrease in volume may not give the same figure with weight reduction. The uniformity feature of the frames is due to the fact that they are constituted from two combined materials: metal and composite.

Based on Tables 5 and 6, it can be observed that practically all frames had their weight significantly reduced (by more than 20%), but only frame numbers 4, 8, and 9 had minor adjustments, with weight reductions of 14, 10, and 17%, respectively. This circumstance is brought on by the frame’s restriction – these frames cannot be optimized further because they have reached ideal size. Frame numbers 8 and 9 cannot be altered significantly (both frames are dominated by a metal structure, while optimization only works for composite material), while frame number 4 is already relatively small; therefore, no significant weight change occurred. Overall, the optimization process has summed up a 26% weight reduction compared to the initial weight.

Table 5

Summary of optimization results on typical frame structure

	Before	After
Frame 1
Frame 2
Frame 3
Frame 4
Frame 5
Frame 7
Frame 8
Frame 10

6.2 Structural strength evaluation

Figure 15 shows the results of the structural evaluation of the overall fuselage construction before the optimization process was carried out. The displacement, plane stress, and shear stress indicate that the fuselage structure shows that the fuselage structure is within the material criteria threshold. The results of the analysis show that the maximum stress values in x-, y-, and xy-planes are 48, 34, and 56 MPa, respectively, which all occur in the area where the metal lugs are installed. In comparison to the metal material strength, the stress is rather low. This evaluation provides the conclusion that the initial structure of the fuselage was designed too optimistically to bear applying loads.

Figure 15

Global FEM results before the optimization.

The topology optimization naturally increases the stress due to the removal of some areas of the frame. This is the expected trade-off which may raise another concern about whether the structure is still safe due to the stress. To evaluate this aspect, global FEM over-optimized fuselage structure was carried out; its results are depicted in Figure 16 showing the structural evaluation after the optimization process. Overall, the value is still under a safe limit. The only concern is the structural response of the metal part in the lug area. Interestingly, the analyses, both before and after optimization, demonstrate that the fuselage frame is not effective in sustaining fuselage load. In reality, the load is dispersed by integral longitudinal structures like girders or stringers. Thus, the maximum stress occurred at the girder in the lug area. In addition, the displacement data shows only 31 mm at the end of the fuselage which is not particularly noteworthy.

Figure 16

Global FEM results after the optimization.

The following findings on the first frame, for instance, are identified by comparing the structural response before and after the optimization. The results show a significant change in shape achieving a weight reduction of up to 42% from the initial weight (Figure 17). Another expected phenomenon is the difference in strength value responses. After the optimization, the stress concentration moves to the removed area. However, the stress value is still below the value of the strength properties of laminated composites.

Figure 17

Stress distribution of fuselage frame before and after the optimization.