Non-Linear Analysis in Post-Buckling Regime of a Tilt Rotor Composite Wing Structure Using Detailed Model and Robust Loading Approach

Abstract

The design and development of a wing for a completely brand-new aircraft represents, in aeronautics, one of the highest challenges from an engineering point of view. The present work describes a novel methodology devoted to execute numeric simulation in a non-linear post-buckling regime to verify the composite wing compliance under the design load conditions. The procedure was developed as part of a wing design and research activity and was motivated by the need to have more realistic results, without standard conservatisms like the no-buckling up to ultimate load, to be of use for achieving further weight savings. To carry this out, it was obviously necessary to ensure that the structural integrity was also guaranteed in the post-buckling regime, above the limit load, and therefore in a highly non-linear regime. The present work illustrates a numerical approach based on non-linear finite element analysis which uses the inertia relief option in order to have a more realistic representation of the structural response of the wing in its real context. All that represents a novelty since, at present, the commercial FE codes allow us to use the inertia relief option exclusively for linear analysis. Obviously, the approach can be applied to any other structural component with similar needs. The obtained results show that the differences between linear and non-linear regime are not negligible and, above all, that it is possible to design a wing (or other structural components) considering, at the same time, the large deformation due to the post-buckling regime, the material non-linearities due to the failures and any other non-linearities in order to achieve the challenging weight requirement of the new aircraft generation.

1. Introduction

Composite structures for aeronautical applications are usually designed to operate under buckling event, beyond the others constraints, resulting in a heavier design the real post-buckled strength is not fully exploited. This is due to the complexity of the stress state that is generated after the buckling onset and the sudden loss of capability to support compressive loads. Furthermore, the non-negligible out-of-plane deformations can lead to overloads in the joint regions with the stiffeners and consequent debonding. Composite laminates, as well, are generally oversized due to transverse shear and normal stresses that may cause early failure due to delamination. In this framework, numerical procedures are gaining more and more importance to investigate and verify the strength of these structures in the post-buckling regime [1,2,3,4,5,6,7], with the final aim to achieve a more efficient structural design [8].

In [9], a design strategy for thin-walled composite structures by considering robustness requirements is presented and traded against numerical results on a composite thin-walled panel by using implicit non-linear finite element analysis solver (Nastran SOL 400, and not SOL 106 due to its limitations. The objective of the above study is to find criteria based on structural energy until collapse, aimed at gaining a weight efficient structure which exhibits a collapse load which is as near as possible to the ultimate load. These criteria, which are currently not part of certification specifications, might lead to an innovative, more robust design with respect to “strength up to ultimate load” criterion and potentially a new approach for the certification of large structures.

For the certification of large-scale structures, a huge number of testing activities is usually performed, at different levels (building block approach) [10], with an overall impact on the time to market and certification costs associated with a new product. Relying on more accurate prediction of the structural behavior of complex structures makes it possible to gain a considerable reduction in both costs and time to market. In this regard, in [7], a methodology based on an appropriate modification of the super-element method, to consider the non-linear behavior, is developed to tackle large scale structure analysis in an efficient way, with a focus on a fuselage section and geometrical non-linearities.

The difficulty in moving from a linear to a non-linear analysis is related to the mode in which loads and constraints are applied to aeronautical structures. In fact, as it is well- known, in the real world, airplanes are not constrained in any way as the aerodynamic forces are balanced by inertial ones. Such an issue is taken into account in the numerical analysis procedures implemented in almost all commercial FE codes but is restricted to linear analysis by using the inertia relief option. Generally, such a restriction does not generate any problems since, as required by current regulations, aeronautical structures must not incur in permanent deformations until the application of the ultimate load and must not present instability phenomena up to the limit load, and therefore there is no real need to simulate non-linear conditions (both material and large displacements). Instead, the case reported in the present work, being focused on a technological demonstrator, can avoid some constraints, preserving the same safety level. In particular, the design rule imposes us to evaluate the stress status by relaxing the no-buckling constraint [11,12], i.e., the instability can be taken into consideration and can occur, but failures must not occur up to the ultimate load. All that, obviously, has the advantage of better exploiting the adopted materials and allows us to achieve further weight savings with respect to the traditional design process. In order to achieve such a goal, it was necessary to extend the normal analysis procedures by including the activation of non-linearities with the inertia relief option, as described below.

The present work is framed within the Horizon 2020 Clean Sky 2 T-WING project, which is aimed at developing, up to testing in relevant environmental (TRL6), the innovative composite wing of the Next Generation Civil Tiltrotor Technology Demonstrator [13]. The methodology used to study the post-buckling structural behavior of the wing is presented, along with main results. Being the vehicle a technology demonstrator and not a final certified product, this study was necessary to assess the need to perform a dedicated sub-component tests and setup a powerful tool devoted to virtual test approach. The phase in which the study was conducted is the Critical Design Review of the project. Before, to reach this maturity stage, other tools devoted to aid the preliminary design of such a wing were developed by the authors in order to achieve the best composite wing structure configuration with the minimum weight [13]. The design of a tilt rotor wing has to respect not only strength requirements but also aeroelastic [14,15,16] crashworthiness ones [17]. All that determined an iterative loop with several disciplines that highlighted that the post-buckling capability was a key factor to achieve the desired weight savings.

The design was mature [18] to invest in heavy calculation by using a Detailed Finite Elements model (DFEM); nevertheless, propaedeutic activity to set up the methodology was already performed during Preliminary Design Phase on coarser models. The Detailed FEM was built for linear analyses (strength and buckling) and when in front of the engineering problem of studying the post-buckling behavior, a number of issues needed to be solved to allow the non-linear calculation to be performed in a reasonable computational time with the available resources.

2. Materials and FE Model Description

2.1. Wing Overview

The wing architecture is the result of the compliance to a number of requirements which are derived from airplanes, helicopters and tiltrotor-tailored airworthiness, and technical and functional specifications. The result is a prescribed position of the spars, ribs and moveable surface geometry, dictated by the Leader Leonardo Helicopters. The wing structural architecture is the result of several optimization processes aimed at assuring strength, buckling and stiffness requirements, with the lowest structural mass [18]. The stiffness requirement, in terms of the first bending and torsional modes, is the most relevant one and it is quite hard to satisfy due to some constraints like the total weight and the accessibility (numerous and wide access holes and removable panels, which tend to reduce torsional stiffness). The wing of the NGCTR-TD has a multi-spar architecture, mainly made by composite solid laminates and sandwich structures, with some other small part in metallic alloys (i.e., the leading edge is made in aluminum alloys). A three-spar wing represents the right compromise in terms of strength and space available inside the wing box. It hosts a number of systems (fuel system, inter-connecting driveshaft, hydraulic and electrical routing) and the interface with the Nacelle Primary Structure at the tips and the fuselage in the central portion.

The detailed Digital Mock-up [19] (DMU) of the wing structure is reported in Figure 1.

2.2. Materials

As previously mentioned, the wing is made almost exclusively of composite materials. Most of the structure is made in fabric CFRP and whose properties are reported. Moreover, some parts have been reinforced with unidirectional CFRP (skin upper panel and stringers) in order to maximize the mechanical characteristics along the principal loading direction. Finally, in order to increase the stability of the ribs, these are made in sandwich with the external skins in fabric CFRP and the core in honeycomb.

The mechanical properties of the adopted materials are reported in

Table 1. Mechanical properties of fabric CFRP-CYTEC 977-2-42-3KT300D-8H-372-1524.

Elastic Moduli–CFRP Fabric
Longitudinal Elastic Modulus-E1 (psi)	811,000
Transversal Elastic Modulus-E2 (psi)	791,000
In-plane Poisson ratio-ν12 (/)	0.05
In plane shear Modulus-G12 (psi)	570,000
Allowable–CFRP Fabric
Longitudinal Tensile Strength–F1T (psi)	113,260
Longitudinal Compressive Strength-F1C (psi)	98,100
Transversal Tensile Strength–F2T (psi)	102,070
Transversal Compressive Strength–F2C (psi)	97,700
In plane Shear Strength–S (psi)	13,233

,

Table 2. Mechanical properties of UD CFRP.

Stiffness Moduli–CFRP UD
Longitudinal Elastic Modulus-E1 (psi)	2,170,000
Transversal Elastic Modulus-E2 (psi)	1,230,000
In-plane Poisson ratio-ν12 (/)	0.3
In plane shear Modulus-G12 (psi)	660,000
Allowable–CFRP UD
Longitudinal Tensile Strength–F1T (psi)	306,000
Longitudinal Compressive Strength-F1C (psi)	194,000
Transversal Tensile Strength–F2T (psi)	6000
Transversal Compressive Strength–F2C (psi)	5100
In plane Shear Strength e–S (psi)	13,500

and

Table 3. Mechanical properties of honeycomb in sandwich structure.

Stiffness Moduli-HC		Allowable-HC
E1 (psi)	1	Tensile 0° (psi)	0.02166
E2 (psi)	0	Compressive 0° (psi)	0.01926
ν12 (/)	0	Tensile 90° (psi)	0.01936
G1Z (psi)	6500	Compressive 90° (psi)	0.06219
G2Z (psi)	3400	Shear in plane (psi)	0.03482

. The reported material properties derived by experimental test campaign, performed by Magnaghi, and they represent the material data adopted for the design; therefore, they include all standard knockdown factors for aeronautical application (statistic, environmental conditions, damage, etc.) [20,21,22].

Moreover, there are several parts, like nose rib, the leading edge, and the tip and root ribs, that are made in aluminum alloy Al7050-T7451 (

Table 4. Mechanical properties of Al7050-T7451.

Mechanical Properties
Young Modulus (ksi)	10,400
Poisson ratio-ν12 (/)	0.3
Tensile strength-ultimate (psi)	76,000
Yield stress (psi)	68,000
Elongation at failure (%)	11

).

2.3. FE Model Description

The Detailed FE model (DFEM), defined in MSC Nastran environment [22], was primarily built to allow for linear strength and buckling verification. It is a global model which takes into account not only the entire wing, but also all the other vehicle sub-components. In particular the fuselage and tails are modelled by means of super-elements (mass and stiffness). The Inter-connecting Drive Shaft, the fuel and the fuel system (including the tanks) are modelled by means of concentrated masses. Being in this work the focus on the wing structure, the control surfaces are modelled in a lumped way, nevertheless their loading effect on the wing is taken into account. The model is built in order to set the equilibrium between the specified air and ground loads with inertia forces, considering each item of mass in the rotorcraft. This is obtained by activating the so-called “Inertia Relief” option and formally constraining the model with a SUPORT in the vehicle center of gravity. Regarding the wing structure, about 95% is made up of shell elements that represent the wing structure and the nacelle primary structure at the wing tips.

The mass of the fuel is distributed in the affected areas of the wing by means of interpolating elements (RBE3), according to the specific load condition.

The entire DFEM model is shown in Figure 2.

In detail, the wing FE model is characterized by the following:

• The 2-D elements (TRIA3 and QUAD4) that are used to model any component (upper and lower skin, leading and trailing edges components, stringers webs and caps, spars, ribs, and moveable surface components). All the parts of the structure have been modelled according to their shape coming from the Digital Mock Up. Aerodynamic surfaces have been modelled from the Outer Mold Line (OML) with an offset to take into account the panel thickness, whereas the internal parts are modelled starting from the middle surface with zero offsets. Figure 3 shows the details of the wing box FE Model.

BEAM elements have been used to model bolted connections [19]. In particular, CBUSH element has been used to model fasteners in the global finite element model. Figure 4 shows the modelling technique used to represent the bolted joint, where CBUSH elements are placed between two RBE2 that are connected to grid points on co-planar meshes.

The CBUSH elements have in-plane stiffness which provides a load path from component A to the component B. Using this modelling technique, the forces from CBUSH elements are quantified as bearing forces in the laminate.

• BEAM, BAR and ROD elements were used to model the kinematics of the moveable surfaces (Figure 5).

  Figure 5. Detail of modelling of kinematics of the moveable surface.

  Figure 5. Detail of modelling of kinematics of the moveable surface.

• RBE3 interpolation elements have the function to transfer loads to the structure (inertial; and aerodynamics).
• RBAR/RBE2 rigid elements to model some links between components (for example the pin-hole connections).
• CONM2 elements have been used to model fuel mass, the structural mass of the different components (material mass density is null), the masses of the cables, etc.

For composite elements properties, PCOMP cards and 2D orthotropic material MAT8 are used whereas PSHELL and MAT1 cards are considered for isotropic material parts. The finite element model can be summarized in

Table 5. Characteristic of finite element model.

Entity Type	Element Description	Number
CBAR	Bar elements	10,781
CBEAM	Beam elements	426
CELAS1	Scalar spring elements	16
CELAS2	Scalar spring elements	126
CONM2	Concentrated mass elements	156
CQUAD4	Quadrilateral plate elements	1,232,122
CROD	Rod elements	20
GRID	Nodes	1,275,487
PLOTEL	Visualization elements	86
RBAR	Rigid bar elements	2
RBE2	Rigid body elements	292
RBE3	Interpolation elements	166

:

The following figures depict the wing model in both composite (Figure 6) and metallic parts (Figure 7). Regarding the composite parts, there are three main difference among them: fabric, unidirectional and sandwich with fabric faces. The metallic parts are all made in aluminum 7050-T7451.

2.4. Load Condition and Boundary Conditions

The total number of sizing load conditions is 50, which includes aerodynamic loads, propulsion forces at nacelle prop-rotor location and inertia loads. They are given in terms of accelerations and rates at vehicle CoG plus the corresponding aerodynamic pressures, balancing loads and concentrated loads such as the prop-rotor thrusts. No boundary conditions, in terms of model constraints, are considered; in fact, since the inertia loads are evaluated via the inertia relief, the model is considered free–free. The perfectly balanced loads are applied to the free–free model and the resulting accelerations, computed by Nastran, are automatically applied to the model generating the subsequent inertia loads, applied to the model in correspondence of the nodes where the concentrated masses are connected to. A representation of the whole set of concentrated loads applied to the model is shown in Figure 8 where the red arrows represent the concentrated loads.

Two sets of load-cases are provided: twenty-five sizing load conditions at International Standard Atmosphere (ISA) ambient conditions plus an additional twenty-five sizing conditions at −45° OAT (Outside Air Temperature) ambient conditions (cold day). The load conditions are used to load the wing FE model and the remaining vehicle in order to perform the inertial relief analysis on the whole aircraft and to take into account the fully balanced vehicle conditions. Each load conditions set includes the most demanding conditions on ground (e.g., taxi, turning, brake roll, and landing), the flight maneuvers both in helicopter and aircraft modes (e.g., symmetric pull up, symmetric push-over, yaw, rolling) and gust conditions. A preliminary analysis highlighted the most critical maneuver for the composite parts of the wing, the load condition so-called LC09 (low-speed pull-up maneuvers in helicopter mode). For the present paper, this loading condition has been proven, but this method can be applied to all loading cases.

The decision to load in a balanced way the entire vehicle is based on a previous experience of the authors and Leonardo partners in order to consider the tails’ balancing loads and the fuselage stiffness at the attachment points. In this, it is possible to reduce the conservativism derived from the different approach to consider the wing as a stand-alone structure constrained at ground at the fuselage attachment points and loaded with aerodynamic, rotor and inertial forces due to the masses (structural and non-structural) inside the wing only. If, for linear static analysis, this approach is very useful, contrary to the non-linear regime, it determines several problems. Basically, this limitation encompasses the entire loading approach. It must be completely revised on one side but it needs to keep the beneficial effects due to the use of balanced loads by means of the inertial relief approach. A method to pass from linear static full-vehicle wing DFEM loaded with balanced loaded to a non-linear full vehicle wing DFEM loaded with equivalent forces is one of the key parts of this work.

In order to obtain, as an output of the inertia relief, the correct aircraft CoG accelerations, the masses of all the components of the model have to be correctly modelled and connected to the structure. For what concerns the aircraft components out of the topic of the T-WING project (i.e., fuselage, tail and nacelles), the WAL has provided global dynamic properties of each component (mass, CoG coordinates and moments of inertia) for each sizing load-case. These data have been used to generate a load-case dependent concentrated mass, representative of the inertial properties of each component, connected to the structure by the use of RBE3.

On the other hand, the mass of the wing has been modelled by the use of several sets of distributed and concentrated masses (structural mass, fuel mass, hydraulic mass, cables mass and fuel systems mass), connected to the proper areas of the structure by the use of RBE3.

In order correctly simulate the inertial contribution of the fuel taking into account the fact that fuel mass acts on different areas of the fuel bladders depending on the local acceleration of the bladder itself, load-case dependent RBE3s have been adopted to connect the fuel mass to the structure depending on the local acceleration of each bay.

For the same reason, as well as the fuel masses, even the masses representing the tank sump plates have been connected to the structure by the use of acceleration-dependent RBE3s. In fact, the fuel storage system design foresees the use of no. 4 retainers for each tank to avoid, under negative g conditions, the collapse of the fuel probe. According to the systems distributed on the sump plates, the lumped masses representing the fuel systems (sump plate, access door and fuel gauging equipment from another research consortium), for each bay, have been attached by means of acceleration-dependent RBE3s to act differently depending on the direction of the acceleration.

The external loads provided by the Leonardo Helicopters have been applied to the model by the use of RBE3s, for each load-case. Each load-case consists of wing aerodynamic loads, nacelle aerodynamic loads, rotor loads and fuselage and tail balancing loads. The wing aerodynamic loads are transferred to wing structure by the use of interface loading grids connected to all the border grids of the upper and lower wing panels, nacelle aerodynamic loads and rotor loads are transferred to the NPS (Nacelle Primary Structure) by the use of RBE3s, and fuselage and tail balancing loads are applied directly to the fuselage DMIG nodes.

As said, the external loads provided by LHD are applied to the model as concentrated forces and moments and the accelerations required to have a self-balanced load condition through inertia loads are calculated by Nastran by the use of inertia relief. The Nastran SUPORT card is defined in the real CoG of each load-case and on all six DoF.

The resulting accelerations, generated by the solution sequence to balance the external loads, are then compared, for each load-case, with the one provided by LHD to check the correct implementation of the whole process.

3. Methodology

Being the finite element model (input of the present activity) best suited for linear static analysis, some modifications have been necessary to allow us to use it in non-linear buckling analysis. These preliminary changes concern the following:

• The correct definition of the offsets for composite elements in non-linear or linear buckling analysis (ZOFFS in 2D elements entries in place of Z0 in PCOMP).
• The original bolts modelling has been simplified by using BEAM elements with appropriate properties.

Furthermore, two main issues have been identified:

• A non-linear analysis does not allow the use of the inertia relief procedure if large displacement non-linearity is activated. A post-buckling analysis requires this type of non-linearity.
• Huge RBE3 elements (close to 100,000 independent grid points in some cases) are used to distribute aerodynamic loads and masses (structure, fuel, …). Note that no density mass is considered for the material used for the several components of the wing.

3.1. Inertia Relief in Post-Buckling Analysis

The structure to be analyzed is not constrained so the inertia relief procedure must be activated in linear static analysis. The inertial loads that balance the applied loads are calculated and added to the load condition; then, the structure is constrained at a point with stiffness in all six degrees of freedom (selected by the user or automatically by MSC Nastran). In these conditions, the static analysis can be performed obtaining the correct stress and strain distribution (depending on relative displacements), and a deformed configuration that depends on the chosen point.

A different mass distribution has been defined in each of the loading conditions to be considered in the analysis, so the SUPORT entry changes for each of them: this point refers to the gravity center of the structure in the specific case.

In non-linear implicit analysis with MSC Nastran, inertia relief can be used in SOL 400, but this approach is only valid for models whose geometry does not change during loading, that is, only in analyses where small displacements are considered (linear analysis).

Note that the old implicit non-linear solution of MSC Nastran, SOL 106, does not sup-port this functionality; this limitation, in addition to others related to, for example, the formulation of elements, the definition of different material models, and the lack of an appropriate simulation of the contact condition led to the decision not to consider it. Considering that a post-buckling analysis needs large displacement/rotation effects, a special procedure has been used to extend the inertia relief functionality in a SOL 400. The basic idea is to use the balanced loads and then constrain the structure in the supported node.

Performing a preliminary linear analysis, it is possible to extract the balanced loads calculated by the inertia relief procedure. To export these equivalent loads, two possible solutions can be considered:

• The creation of an appropriate DMAP procedure that uses a module able to export the data available in a load vector directly in the format of the FORCE and MO-MENT entries.
• Execution of the linear static analysis with inertia relief in whose input file the request for the output of the applied loads (OLOAD = ALL in the Case Control Deck) and for storing the data in a file in HDF5 format (MDLPRM, HDF5, 1 in the Bulk Data Deck) are present. An appropriate external code allows accessing to the file in HDF5 format, creating the resulting FORCE and MOMENT entries.

It was decided to use the second approach. In fact, in the first case, the FORCE/MOMENT cards are written in the large field format (16 columns for each field) and the numerical values are represented in exponential form. In this representation, four characters are engaged for the exponential term, thus effectively reducing the number of significant decimal digits with consequent possible loss of precision.

Accessing to the HDF5 file through an external code created ad hoc, it is possible to manage the format of the numerical value to insert in the appropriate field of the FORCE and MOMENT entries; this means that it is possible to use a floating point representation of the number allowing us to consider more significant digits in the 16 characters. The Figure 9 shows the flowchart to define the equivalent inertia load to be applied.

Comparative tests have been performed between the results obtained in the linear analysis that use the standard inertia relief procedure and the one that considers the structure constrained in the SUPORT node and the equivalent calculated loads applied to the proper grid points. This test was successful and, therefore, this load can be safely used in the subsequent analysis without the need to implement the inertia relief procedure. The generated FORCE/MOMENT entries have been used as the load condition in a non-linear analysis in which the node, originally declared in the SUPORT entry, has been constrained in all its degrees of freedom. In this way, we have reproduced what is performed in the linear static analysis with inertia relief but having the possibility to set a step-by-step solution in the SOL 400 with large displacement non-linearity activated (PARAM, LGDISP, 1). In this condition, the post-buckling analysis can be performed.

3.2. Issue Related to Huge RBE3 Elements for Non-Linear Analysis

The finite element model of the wing RBE3 elements have been used to distribute loads and masses. This choice has led to the definition of many RBE3 elements characterized in some cases by an extremely high number of independent nodes; this determines in certain types of analysis a considerable loss of efficiency.

In fact, in modal analysis, the high density of the mass matrix leads to an excessive increment in the resources required for the eigenvalue extraction (in terms of space and time). In the same way, in non-linear analysis, the Lagrangian representation of the de-pendency relationships leads to consider the RBE3 as a real structural element that, in the presence of a high number of independent nodes, determines a considerable increment of stiffness matrix density with a consequent reduction in the efficiency of the calculation. Even in case of a linear formulation for RBE3, there is an increase in the time spent in some DMAP modules (MCE1 and MCE2); this does not introduce performance problems in a linear static analysis as these operations are performed only once. In contrast, in a non-linear analysis, these modules are recursively performed at each iteration resulting in a significant increase in calculation time.

Some tests executed on the finite element model of the wing containing these RBE3 elements confirm the fact that this modelling solution is not suited for a non-linear analysis. In fact, for one of the loading conditions, the solution was obtained after a little bit more than 72 h (3 days), while in another one, the calculation was stopped after 20 days during which the convergence was reached only up to the 80% of the total load. The very high reduction of the performance in this last load condition is due to the presence of a higher number of huge RBE3 (for fuel mass distribution), which are not present in the first one. Both the compared analyses were run on the same Linux server (256 Gb of physical memory and 28 CPUs) with the same launch settings (memory requests, number of CPUs, …).

The solution has been found in the use of the MPCFORCE output. In fact, in this way, it is possible to evaluate how the equivalent load applied to the dependent node is distributed in the independent nodes belonging to the single RBE3. Since many nodes may be declared as independent in more than one RBE3, it is not possible to extract this information in one shot.

In the definition of the model, the RBE3 elements have been distributed in different files and grouped in order that each file contains all the rigid elements that describe a specific feature. For example,

• A file contains all the RBE3s that describe the aerodynamic load;
• Another file contains the masses of the structure (remember that no mass density is defined for the structural elements);
• A third file contains the masses of the fuel.

The task was therefore performed by considering all the large RBE3 elements contained within each of these files. Starting from this idea, the following procedure has been built (Figure 10):

• Considering the complete wing model, a preliminary inertia relief job is performed to create an HDF5 file that contains the output request for the balanced loads (OLOAD = ALL and MDLPRM, HDF5, 1).
• A specific Python code reads the HDF5 file and create a file in which the FORCE/MOMENT entries representing the balanced loads are written.
• The complete finite element model of the wing is imported and a group containing only the RBE3 elements defined in each of the mentioned files is created.
• Each of these groups is populated with the structural elements (and the connecting grid points) that are attached to the independent grid points of the specific set of RBE3.
• The independent grid points of the RBE3 elements belonging to each of the created groups are constrained.
• Each of the dependent grid point of the RBE3s is loaded by considering the corresponding balanced load (aeroelastic or inertial) obtained in the preliminary inertia relief analysis.
• For each group, a specific MSC Nastran input file is exported. The output request for MPCFORCE to be stored in the HDF5 format is activated
• For each of the generated input files, an MSC Nastran job is run, creating several HDF5 files containing the MPCFORCE output associated with each set of RBE3 elements.
• Read each of these HDF5 files and extract the MPCFORCEs in the format of the FORCE and MOMENT cards representing the loads distributed in the independent grid points of each set of RBE3 elements by using an external properly created Python code. Furthermore, considering other non-huge RBE3 elements distributing loads and masses, the file containing the balanced loads in the FORCE/MOMENT format is modified by eliminating the entries that refer to the dependent grid points of the huge rigid elements.

By using this procedure, it is possible to eliminate the large RBE3 elements by re-placing them with the equivalent loads reported at the independent nodes. It involves the redefinition of the files containing these RBE3 (eliminated), for example copying them in new files renamed by adding the suffix ‘No-RBE3’ to the original name.

3.3. Semi-Automatic Procedure to Prepare Data for Post-Buckling Analysis

According to what reported in the previous section, a semi-automatic procedure was set up to prepare the data necessary for the post-buckling analysis of the wing under different loading conditions. The procedure is summarized in Figure 11.

Each of the two external codes require a specific input file containing the following information:

STEP 1 Procedure

• The name of the HDF5 file containing balanced loads data.
• The name of the file to store the FORCE/MOMENT entries equivalent to the inertia relief loads.
• The identifier for the above FORCE/MOMENT entries.
• The name of the Patran session file containing all the commands needed to build the dummy models and create the MSC Nastran input file for each of them.
• For each set of huge RBE3 elements.

  ○

  The name of the file containing the current set of RBE3 elements.

  ○

  The name of the group associated with the current set of RBE3 elements.

  ○

  The name of the file containing the FORCE/MOMENT entries defining the balanced loads associated with the current set of RBE3 elements.

  ○

  The identifier for these FORCE/MOMENT entries.

  ○

  The name of the files in which store the FORCE/MOMENT entries representing the loads distributed in the independent nodes of the huge RBE3 (one for each huge RBE3 file).

STEP 2 Procedure

• The number of files with huge RBE3 elements.
• For each set of huge RBE3 elements.

  ○

  The name of the file that contains the list of the dependent grid points of the huge RBE3 elements (It is generated in the STEP 1 and has predefined name).

  ○

  The name of the HDF5 file generated by running the input file relative to this set of RBE3.

  ○

  The name of the output file generated by the code containing the FORCE/MOMENT entries that describe the loads distributed to the in-dependent grid point of this set of RBE3.

  ○

  The identifier for these FORCE/MOMENT entries.

To simplify the generation of all, a special Graphic User Interface (GUI) has been created (Figure 12). This was designed based on the organization of the structure data and an output format obtainable from a calculation performed with MSC Nastran; as it is, the procedure and the interface in particular are therefore usable within this specific project, although the general idea could be applied to any other situation since it is all conducted in Python.

Using this GUI, it is possible, for example, to access the folders in the disk to select the files or to use some predefined defaults for some input data. Note that, for this project, the several input files with ‘congruent’ behaviours have been organized in different folders with assigned and unmodifiable names. Considering that the huge files are all in specific folders, the procedure allows us to choose between them to directly access the proper location. The output of this GUI consists of three files:

• The input file for the Python code related to STEP 1.
• The input file for the Python code related to STEP 2.
• The batch file containing all the commands that allow you to perform the following operations in sequence:

  ○

  Run Python code for STEP 1.

  ○

  Execute in batch in a MSC Patran session to create the input files for the dummy models associated with each of the set of huge RBE3 elements.

  ○

  Execute a MSC Nastran job for each of the generated input files.

  ○

  Run Python code for STEP 2.

At the end of the execution of this batch script, several files, containing the loads distributed at the independent node of the huge RBE3 elements, are stored in a specific location. These files are ready to be included in the post-buckling input file by using the INCLUDE command.

3.4. Definition of the Contact Conditions

The provided finite element model shows, locally, some shapes not compatible with the real situation due to penetration or gaps between adjacent components. The bolts as-sure the connection between the components but allow for relative displacement in both directions (traction and compression); this means that the two components can move away (allowed displacement) but can also approach each other (condition not allowed if in contact).

Likewise, in buckling analysis, the possibility of relative motion in the normal direction between the two bodies connected by bolts leads to the identification of unrealistic instability modes. In fact, the only possible modes are the ones in which the two bodies move together. The presence of these unrealistic modes introduces some related problems:

• Difficulties to identify unrealistic modes;
• Need to compute useless modes More modes mean more computation time.

This leads us to model the contact conditions to better simulate the real behavior of the structure. Depending on the analysis type, different approaches have been considered:

• Static Analysis

Considering the extension of the connection between stiffeners and skins (lower and upper), it has been decided to use a glued contact condition for them. Conservatively, the bolts assure the connection between the other adjacent components.

• Buckling Analysis

In this type of analysis, a preliminary static analysis must be performed to define the differential stiffness associated with the specific load condition. Then, this matrix is considered in the eigenvalue problem that allows calculating the instability modes. In MSC Nastran, this means that two consecutive SUBCASEs must be considered.

In accordance with the static analysis approach, an inertia relief problem is set and a glued condition is considered for the Stringers–Skin connection. Then, in the buckling analysis SUBCASE, an MPC is added to consider the glued contact conditions of all the other adjacent components or only those ones that introduce unrealistic instability modes in the range of interest.

The MPC relationships related to glued contact conditions are associated with all or part of the adjacent components by performing a check run (PARAM, CHECKOUT, YES) in which the command ‘NLOPRM MPCPCH = BEGN’ is specified in the Case Control Deck.

Then, the SOL 105 Case Control Section is defined as shown in Figure 13:

The presence of the inertia relief in the static analysis leads us to introduce, in the buckling analysis, the constraint in the SUPORT node for congruence with the condition in which the differential stiffness matrix is calculated.

• Non-linear Analysis (Post-Buckling)

In this model, about 150 contact pairs have been defined. Considering a touching contact condition for all the contact bodies in a finite element model like the one used in this activity could lead to large resources requests in terms of computational time and disk space. In the initial job all the contact bodies have been considered in glued contact condition. Then, depending on the structural response of the different parts of the model and on the buckling modes that fall within the range of the applied load, for certain contact torques, we passed to the touching contact condition.

Contact analysis requires the correct (or the best possible) definition of the geometries of the connected components; in fact, in non-optimal conditions that overestimated or under-estimated the results may occur because the relationships introduced do not correctly represent the real behaviors of the structure in these regions. When defining the contact conditions between components described by 2D elements, both the thickness (the contact surfaces are created at the top and the bottom of the element) and the offset (repositioning of the middle plane of the element with respect to the geometric position of the nodes that define it) are involved.

The last FEM received solved some, but not all, issues related to PCOMP and the offset definition. Penetrations or gaps still exists because the finite element model was initially created using a geometry that in some parts is no longer compatible with its current setting.

Due to some geometrical misalignments, internal constraints were generated. For such a reason, it was not possible to use the clearance functionality. Figure 14 shows some details on the geometrical misalignments with the rendering of the shell thicknesses.

The contact tolerance parameter ERROR has been used to resolve these geometrical issues. It defines the region across the surface contact in which the contact condition is detected.

Refer to the standard MSC Nastran documentation [21,22] for details and explanations related to this specific topic.

To achieve this contact setting, some modal analyses on the free–free structure with various combinations of contact tolerance values were performed. The analysis of the contact status provided the tool to modify the parameter in those areas where adjusting the ERROR parameter unrealistic modes could be eliminated. Then, the frequency value of the rigid modes has been compared to the first flexible mode to confirmed that the contact conditions is not introducing internal constraints.

As shown in

Table 6. Distance tolerance (ERROR) definition in the contact pairs.

Distance Tolerance (Error)	% of the Total Contact Bodies
0.05	49.5%
0.1	34.0%
0.15	10.7%
0.2	5.8%

, approximately 85% of the contact pairs require a distance tolerance of 0.1 and therefore a value much lower than the original setting of 0.5 used in the first FEMs. In any case, the positive effect of this setting is that, for any contact pair, we have avoided the use of a very high value for ERROR., even if the best approach should be to not use it. Such data could help us understand where to work to improve the geometric representation in the finite element model. It has been considered acceptable, in the mean-time, that an updated model will be available.

3.5. Miscellaneous

Other changes to the original model and input file have been made to overcome both existing limitations in a certain type of one-dimensional elements and to overcome some problems verified in the analysis of the results of some preliminary non-linear analysis.

• BAR elements in non-linear analysis

In the finite element of the wing, BAR elements are extensively used to model some mechanisms. Furthermore, in non-linear analysis, bolts are also modelled using this kind of elements.

In MSC Nastran, native BAR elements do not support large displacement, so it was necessary to transform them to be used in a non-linear analysis. Several actions are possible:

○

Use of large strain elements by adding the PBARN1 entry to the standard property entries (PBAR and PBARL). In this case, it is convenient to also change the other elements. In general, this solution has been successful, but in this case, it has been discarded because the tests performed showed convergence difficulties and therefore a reduction of the computational efficiency.

○

MSC Patran could be helpful in modifying the element type activating the option “Convert CBARs to CBEAMs” in the input file generation phase. It has been also discarded as it requires the replacement of all the cards associated with the BAR elements within the various files used to describe the model. There is a high possibility of errors in this operation.

○

Using the entry “MDLPRM, BRTOBM, 1” in the MSC Nastran input file, the con-version of BAR elements to BEAM is automatically activated. This solution has been used.

• Some modification in connections in the flaperon region

Some non-correct modeming solutions were highlighted in the flaperon area during the execution of the preliminary post-buckling analysis. The affected area is the one shown in Figure 15.

The modifications are relative to the modelling of the hinge or the spherical joint that connects the flaperon shaft to the end supports.

For example, in one of these connection regions, two coincident grid points have been created; one of them is the independent grid point of the RBE2 elements representing the hole of the support while the other one is the end node of the BEAM that defines the pin. Another RBE2 element connects all degrees of freedom of the two coincident nodes minus the rotational one around the hinge axis (in case of spherical joint all the rotational degrees of freedom are not connected) as reported in Figure 16.

4. Results

This section shows the main results concerning the most critical load condition and therefore the sizing one for the entire wing, i.e., the LC9 which, respectively, represents the symmetric pull-up manoeuvre.

In order to show the advantage in using the developed procedure, both the results of the static linear analyses (linear static and eigenvalue buckling) and the results of the non-linear analysis with the application of the methodology are presented.

Figure 17 reports the deformation state at limit load related to LC09. As expected the global deformation is upward and it is symmetric. The maximum displacement is localised at the wing tip and it is equal to 5.97 in.

Table 7. Linear buckling eigenvalues.

LC 09
Mode ID	1	2	3	4	5	6	7	8	9	10
Eigenvalue	1.395	1.499	1.527	1.564	1.703	1.711	1.763	1.782	1.818	1.909

reports the eigenvalue of the linear buckling analysis and Figure 18 shows the eigenvector related to the first two modes whose scaling factors are below the ultimate load value. The buckling affect only the wing tip ribs, in particular, for the first two modes less than 1.5 (that correspond to the ultimate load).

Using the developed non-linear procedure, the global deformation state is very similar to the linear one and the difference is about −3.9%; thus, the maximum displacement at the wing tip is 5.73 in. about the buckling. In this case, it is not possible to identify a well-defined value. Figure 19 reports some frames of the non-linear analysis in order to highlight a load range for the buckling onset. Up to 80% of LL no instabilities can be observed. Beyond such a value, a local buckling that involves the upper skin in the central box and in the trailing edge starts. For sake of clarity, no displacement contour plot has been shown since the global displacement at the wing tip is bigger than the one related to the local buckling at the wing central box after the rear spar. Through the analysis of the deformed shape, it is possible to note the local buckling region.

Linear buckling results for LC09 suggest that two buckling events affect the upper skin of the wing. To better understand this phenomenon in the non-linear regime, local deformation at several load increments is reported in Figure 20 for the upper skin. Results suggest that a local deformation on the trailing edge of the upper skin panels happen. However, such a deformation seems to be very gradual and so it cannot be related to buckling even up to 150% of limit load. This means that the non-linear approach allows us to consider some effects that shift the instability phenomena to levels above those considered safe.

By running a non-linear analysis, it was possible to evaluate the stress state load increment by load increment and to apply a progressive failure approach (less conservative than the first ply failure approach). The failure criteria chosen to identify the onset of the ply failure is the Tsai-Wu. In order to correctly estimate any possible failure, the elements close to the bolt elements have been not considered since, obviously, they represent a non-realistic stress concentration. Figure 21 reports the envelope of the failure index on the entire wing structure and it is possible to note that the maximum value is equal to 0.609; therefore, no failures occur (the envelop represents the maximum value among all plies contained in an element).

Figure 22 shows the von Mises stress distribution on the metallic parts, that, as stated before, are made in aluminium alloy. The maximum value is less than the yield stress and therefore a positive margin of safety can be determined both for yield and for ultimate. There are only a few elements with a bigger stress value but all are related to the bonding and bolted conditions (contact issue) and therefore they can be neglected (such elements have been removed by the contour plot and it is clear by the next figure).

5. Conclusions

This article deals with the development of a procedure for FE non-linear analysis in the post-buckling regime using detailed models (high computational cost). In particular, the methodology was applied for designing the wing of the NGCTR, mostly made of composite material, but it is applicable to any structural components with similar needs. Generally, these types of structures are analyzed in the linear regime since it is not expected to reach large deformations (especially out-of-plane due to local and global instabilities); therefore, the computational cost is still acceptable. In the present case, instead, there is the need to exploit the maximum capability of the composite materials, pushing them to work also in a post-buckling regime and up to ultimate loads. In fact, the no-buckling limit has been moved from the limit load to the ultimate load. Thanks to the non-linear analysis, it was therefore possible to verify that the margins of safety are still positive even after the buckling onset, which occurs between the limit load and ultimate load, and all that can lead to non-negligible further weight savings. The procedure therefore dealt with two distinct problems: the extension of the use of the inertia relief option in a non-linear regime and the management of REB3 with a huge number of independent nodes. The RBE3 issue is linked to the loading procedure adopted especially in the aeronautical field which allows the points of application of inertial and mechanical loads to be managed more effectively. The reported results demonstrate that developed non-linear procedure does not significantly affect the global response of the wing; indeed, the difference between linear and non-linear wing tip displacement is less than 4%, as expected since no local and global failure occurs up to ultimate load. But, by using the non-linear procedure, it was possible to fully study the instability phenomena and it was possible to evaluate, in further detail, the stress status re-distribution after the buckling onset. Furthermore, by using the non-linear approach, it is possible to adopt the progressive failure procedure, and thus the composite material can be exploited for its real potential and maybe all that could be used for further structural optimization.