Aeroelastic Response of Aircraft Wings to External Store Separation Using Flexible Multibody Dynamics

In aviation, using external stores under the wings is a common method of carrying payload or fuel. In some cases, the payload can be rigidly attached to the wing. However, stores must often be ejected during flight for aircraft, such as military type, which carry drop tanks and missiles. This may cause the wing to respond dynamically with increasing amplitudes, due to the impulsive load of ejection and the change of total mass. This is especially critical in aircraft with highly flexible wings, such as those with high aspect ratios. In this case, it is crucial to evaluate the wing response to store separation, which requires a suitable simulation environment that is able to support nonlinear and multidisciplinary analysis. To address such a need, this work presents the use of flexible multibody dynamics in the simulation of wing response to store separation. To demonstrate, a highly compliant wing was selected with a rigid body that was mounted on the wing to represent an external store. The time marching simulation of the wing before and after the store separation was presented to show the features and benefits of the method.


Introduction
An external store is a component that is required when there is no available space inside an aircraft's body. A classic application of external stores is to the weapons or drop tanks used in combat aircraft, which have been very common since the early days of aviation [1]. Recently, unmanned aerial vehicles (UAVs) have also started carrying wing-mounted stores, which are typically weapons for combat purposes [2]. With an increasing future interest in UAVs [3], especially due to them being a low-cost alternative to conventional combat aircraft [4], it is expected that the adoption of this role will continue to increase. In brief, manned and unmanned aircraft both carry external stores for ejection at some point during a mission.
These external stores are usually positioned under the wing, alhough sometimes under or above the fuselage is preferred. Depending on the store weight, size, and geometry, store separation causes two consequences: (i) a change in aircraft configuration and, hence, wing behavior, and (ii) the transient response of the wing and store during the short period of the separation process. Whilst the first case can be considered to be two distinct aircraft, the latter should consider the two configurations. Therefore, the store ejection has an intrinsic characteristic that should be handled cautiously. In particular, the separation of external stores from an aircraft wing during flight is a safety concern. An obvious issue is the strike of the store after being released. To address this, the trajectory of the stroke needs to be carefully estimated. This is usually achieved through conducting high fidelity computational fluid dynamics simulations [5,6]. In most safe separation studies, the wing is assumed to be rigid. Still, there are some works that consider the effect of wing elasticity on the store separation process [7][8][9]. In these studies, the focus is on the store trajectory, rather than the response of the wing.
In addition to the wing-store strike, the response of the wing to the store separation due to impulsive load induced by the ejection and change in the mass distribution carried by the structure is another safety issue. There are no reported incidents that are related to wing response to store separation, but recent works show that there is an increased interest in this problem due to the prevalence of more compliant wing structures. For example, dynamic loads that are associated with local impulsive excitation are obtained using mode displacement and force summation methods in refs. [10,11] investigated the aeroelastic behavior of wings due to store ejection using fluid-structure coupling. The structural model is a linear model that is obtained using a Finite Element Method. Additionally, a frequencydomain method using linear FEM is used in Ref. [12] to compute wing response. Linear structural models could provide accurate results for relatively stiff wings. However, as wing compliance increases, this could lead to moderate and large deformations. As a result, non-linear behavior may arise due to changes in structural, dynamic, and aerodynamic loads. Therefore, the resulting wing response could be very different from that of a linear structure. As in the case of any non-linear dynamical system, the wing could experience limit cycle oscillations or chaotic response [13,14]. Therefore, the assumption of linearity can result in misleading results as the non-linear effects due to the wing flexibility increasing. A better approach would be to use a method that can support nonlinear simulations.
The problem of store separation is, in general, time-dependent and non-linear. Such problems could be handled using spectral methods [15,16] or time marching simulations. This work prefers the latter and it presents the use of flexible multibody dynamics to simulate the wing aeroelastic response to store separation. When considering the efficient handling of large deformations and rotations and the ease of simulating non-linear timedependent problems, flexible multibody dynamics is a promising method in solving aircraft store separation problems. To the author's knowledge, this work is the first utilization of multibody dynamics in the stability assessment of aircraft wings releasing external stores during flight. A solver with an extensive aeroelastic library was selected to achieve realistic results, which includes geometrically exact beam formulations, and static and unsteady aerodynamics. The aeroleastic response of the wing is analyzed using time marching simulations of the wing, both with and without store separation.
The contents of the paper are organized, as follows. Section 2 describes the method and presents a verification of the method using a simple mechanical system. In Section 3, a high aspect ratio aircraft is explained to show the basic elements that are used to model the problem. Section 4 presents the analysis results and demonstrates the use of multibody dynamics in the dynamic response analysis of a high aspect ratio wing. Section 5 concludes the paper.

Method
This section describes the method of simulating a high aspect ratio wing response to store separation and then verifies its use for the separation aspect using a two degree of freedom model.

Multibody Dynamics
Aircraft wings are lightweight structures, and some are even specially made compliant to increase payload carrying capacity or flight efficiency. As a result, wings can demonstrate non-linear behavior with moderate to large deformations. Additionally, store separation changes the configuration of the simulation model, which requires special attention. Therefore, the method that is used to compute the wing response to store separation should be a non-linear and comprehensive solver. This work uses MBDyn open-source multibody dynamics code [17]. MBDyn has been developed at the Aerospace Science and Technology Department of Politecnico Di Milano and it contains a wide range of elements, including an extensive fixed and rotary wing aeronautics library. MBDyn simulates heterogeneous dynamical system that are based on first principles. In other words: differential algebraic equations are directly solved [18]. MBDyn can provide the following capabilities to its users: • nonlinear dynamics of rigid and elastic bodies; • geometrically exact beam and shell elements; • lumped and super elements; • smart materials, electric and hydraulic networks and active control; • essential fixed-wing and rotorcraft aerodynamics; and, • links to external solvers for co-simulation of multiphysics problems, such as Computational Fluid Dynamics (CFD) and terradynamics.
Being an open-source code, MBDyn is extensively used in research and development. Some examples include aircraft maneuver simulations with coupled CFD [19], the simulation of tiltrotor flight mechanics [20], aeroservoelastic analysis of wind turbines [21], and the trajectory optimization of robotics [22].

Verification
Multibody dynamics is a well-established discipline of computational mechanics and it has been verified and validated extensively against problems of distinct characteristics and variable complexity. Still, the separation of two bodies using MBDyn has not been done, and there is a need to verify the approach on a simple problem before applying it to a complex multidisciplinary problem regarding the separation of two bodies. To demonstrate the capability of the method in solving store separation, experiments that are conducted in a wind tunnel would be the most suitable way. In the absence of the necessary resources, an alternative would be to verify against an analytical solution. The latter approach is acceptable, since the simulation tool has been verified in many complex problems [23][24][25] and only the separation aspect of the problem needs further attention.
For this purpose, a two degree of freedom system mass-spring-damper lumped parameter model was selected, as shown in Figure 1. In this form, the two degree of freedom model has an analytical solution and, hence, can be easily compared to a solution obtained using multibody dynamics. Figure 1 presents the schematic of the verification problem. The problem consists of two masses, one of which (m 1 ) is connected to the ground by a spring (k 1 ) and damper (c 1 ). The grounded mass is connected to another mass (m 2 ), which is free on its other side. In the initial condition, these two masses were interconnected by spring (k 2 ) and damper (c 2 ) elements. There is harmonic forcing on (m 1 ) at a 20 Hz excitation frequency. At 0.5 s, the interconnecting spring and damper connection was deactivated to observe the response before and after the separation. Equation (4) gives the equations of motion, with the numerical values of the parameters being reported in Table 1.

Parameter Value
Mass (m 1 , m 2 ) 7000, 500 (kg) Stiffness (k 1 , k 2 ) 1 × 10 8 , 1 × 10 8 (N/m) Damping (c 1 , c 2 ) 1 × 10 5 , 1 × 10 5 (kg × rad/s) Force on main body 3750 (N) at 20 Hz. Time of store separation 0.5 s The lumped parameter system shown in Figure 1 was modeled in the MBDyn environment to demonstrate how a similar analysis can be conducted in a multibody software, proving its validity by comparing the simulation results with that of the analytical solution.
In the MBDyn framework, this system consisted of two rigid structures, three structural nodes, five joints, and one force. The 'total pin' joints were located between the ground and main body and store and ground, to constrain the two translational and three rotational motions. Two 'rod' joints were used to represent the linear spring and dampers. An 'absolute force' was applied in the form of a sine wave at 20 Hz on the main body. In the store separation analysis, the mass was released at the specified time. The 'driven' element was used to separate the store from the main body at a given time that is prescribed in the input file. The multibody model was simulated using MBDyn, and the resulting displacements, velocities, and accelerations were recorded. Figure 2 shows the simulation results for both rigid bodies. Before the separation, they both have almost identical motion due to the very rigid spring connection. After the two masses were separated, the forces acting on the released mass (m 2 ) vanishes, and it follows a constant velocity motion, which is exactly equal to the velocity at the instant of separation. On the other hand, although there is still forcing applied on the (m 1 ), the amplitude of its motion changes. The resulting displacement amplitudes of (m 1 ), which are obtained independently while using the analytical solution and multibody dynamics simulation, are compared in Table 2. The reported values are identical, which verifies the discontinuity aspect of the separation problem.

Demonstration Model
This section describes the aircraft wing that was used in the analysis. Firstly, the aircraft is described, and then the multibody dynamics realization of the aircraft wing model is explained.

Aircraft
The proposed approach can be used in any type of aircraft and wing. Having said this, highly compliant structures can benefit from the method more than relatively stiff wings, as the former types are likely to exhibit nonlinear behavior due to moderate and large deformations. For this purpose, a high aspect ratio was judged to be more suitable as an example, since it is expected to be highly compliant to sudden time dependent forces, as experienced by the wing after the release of an external mass. As a result, the benefit of using a fully nonlinear solver could be more easily demonstrated. Consequently, an unmanned aircraft with a high aspect-ratio wing was selected as the demonstration vehicle, as presented in Figure 3. The aircraft of interest is a high-altitude long-endurance (HALE) unmanned aircraft and it has been previously studied for generalized eigen-analysis in Ref. [26] and for limit cycle oscillations in Ref. [13]. Having highly compliant wings, this aircraft is a good candidate to demonstrate the wing response to store separation. Because this work focuses on the wing response, the fuselage was not included and a single cantilever wing was considered, with the numerical values that are reported in Table 3. An external store under the wing is considered with the parameters that are given in Table 4. The wing is simulated at a high altitude flight condition, the conditions of which are available in Table 5.  Table 4. Numerical values of the store.

Parameter Value
Mass of Store 10 kg Store Support Stiffness 1 × 10 8 N m −1 Store location with respect to the wing root 7 m

Parameter Value
Altitude 20 km Air Density 0.089 kg m −3 Flight Speed 40 kn

MBDyn Model Description
The wing model without the store, and its separation process, was previously developed while using MBDyn [17,26] with the elements that are presented in Figure 4. Additional features were added to this model to simulate the wing response to store separation. Namely: stores, the store-wing connection, and drive element to release the store. The final overall multibody dynamic model consists of the following:
Structural model: 16 three-node finite-volume beam elements [27] defined the 16 m span wing structure. The whole wing structure required 32 nodes for structural discretization, each was assigned inertial properties.

3.
Beam cross-section: properties were given as a diagonal (6 × 6) constitutive law matrix, which remains constant along the span. 4.
Aerodynamic model: the wing was discretized by aerodynamic elements formulated by lifting line theory. The aerodynamic elements included both steady and unsteady aerodynamics, therefore helping to achieve realistic static deflection and dynamic response.

5.
Boundary condition: ground joint was used to support the wing root. 6.
Store: a rigid body at X store = 7 m spanwise location was defined to represent the store, which was connected to the wing via spring and damper, as shown in Figure 5.
The store has no aerodynamic model attached to it. 7.
Separation: the driven element was used to keep the store connected to the wing for t < 25 s. 25 s was defined based on the time required to achieve the steady-state solution of the same wing under the same flight conditions. 8.
Gravity: introduced to achieve realistic wing static response of a highly flexible wing and motion of the store after separation.  It should be noted that MBDyn has an extensive library of mechanical parts, joints, and forcing options [17]. Additional elements could be added to make the model more complex. For example, the rod connection could include dedicated joints, friction, or an impulsive force could be included at separation. However, within the scope of this work and the available data, the above multibody dynamics models can be considered to be sufficient in terms of the essential components.

Results and Discussion
This section presents the results that were obtained using the model described in the previous section to demonstrate the benefits of using multibody dynamics in aircraft wing response to store separation.

Simulation without Store Separation
The first results present the simulation of the wing-store configuration to explain the distinct characteristics of the simulation. The displacement of the wing tip was obtained from the simulation and is shown in Figure 6. The simulation starts at t = 0 s and the tip starts from a reference position. When the simulation starts, the gravity and the flight speed are introduced gradually, rather than instantly, to aid with convergence. Gravity dominates just after initiation as the wing flaps down, being mainly resisted by the elastic forces. However, as the aerodynamic forces build up with increased airspeed, the wing recovers and the wingtip moves upwards. Gravity, aerodynamic, and elastic forces come to equilibrium near t = 20 s, after which the wingtip position remains constant.  Figure 7 presents the response of the wingtip before and after the store separation, which is compared to that of the 'no separation' case. At t = 25 s, the store is released from the wing. As a result of weight change, an impact load acts on the wing and the wing tip moves up and starts oscillating. The oscillations are unstable in nature with increasing amplitude, while the 'no separation' case remains constant, which implies a stable blade. This result could be attributed to the non-linearities of the whole system as a result of the change in wing state after separation. It is also worth looking at the deformed shape of the wing. Figure 8 shows the displacement at the beam nodes, which, in turn, defines a displacement shape. Initially, the wing nodes form a straight line at t = 0 s. As the gravity dominates, but is resisted by the elastic forces, the wing flaps down and comes to a lower state around t = 5 s. Subsequently, the aerodynamic loads build up and, consequently, the wing flaps up. This continues until the wing reaches steady-state. Right after the separation, the wing starts an instant upward movement. As the time increases, the amplitude of wing displacement gains amplitude and the solution diverges eventually. With Figure 8, it can be observed that, as the amplitude increases, the oscillatory component of the motion is not visible. In order to better observe the time-dependent motion that the wing undergoes, the oscillatory wing deformation was obtained by subtracting the steady wing motion just before the separation from the total displacement. Consequently, Figure 9 presents wing oscillatory motion after the store is released. The displacement shapes were identified at different time steps within the same cycle (near 50 s). The curves suggest that the oscillatory motion is a combination of the typical first and second bending modes of a cantilever beam.

Stability Enhancement
Two analyses were conducted in order to demonstrate how the suggested method can be used to improve the stability of the wing motion after store separation.

Wing Stiffening
Increasing structural stiffness is an obvious choice for stabilizing the wing. To achieve this, a symmetric rectangular thin-walled beam was assumed inside the wing, as shown in Figure 10. This could be seen as an oversimplified assumption, since aircraft wings are much more complex than a box beam representation. Nevertheless, considering that the beam constitutive matrix was symmetric and constant throughout the wing in the original model, a simple box beam assumption could be acceptable within the scope of this work. It should also be noted that the proposed method could be used in parallel to advanced wing geometries as long as the cross-section constitutive matrix is provided. However, the objective here is to estimate the required amount of increase in the stiffness of the wing, leading to reductions in the oscillations after store separation. For this reason, the constitutive matrix needs to be modified. Consequently, the change in the box-beam thickness was related to the change in the area, moment of the cross-section area, and wing mass.  Figure 11 shows the tip response of the wing in its original stiffness and modified version after increasing the thickness. Whilst the upper plot presents the full simulation, the lower plot focuses on the oscillatory response by focusing on the final time steps. It is clear from the lower plot that the oscillations of the wing reduce as the wing becomes stiffer with an increasing thickness. It can also be observed that there are no oscillations when the box-beam thickness of the wing is increased by 6%. Such analyses, with more sophisticated finite element models of the cross-section, can be helpful, when a wing is not initially designed for store separation, but later decided to do so. The amount of wing design changes that are required to preserve the stability characteristics could be identified.

Vibration Control
While structural modifications can enhance stability, the wing can also be equipped with sensors and actuators to control the vibrations [28]. MBDyn can support such an analysis through its control elements. There is no linearization requirement and the controller can take effect during the time marching simulations. For this purpose, the proportional-integral-derivative (PID) control was implemented. Based on the measured acceleration at the wing tip, a control moment was applied to twist the blade and attenuate the oscillations, as shown in Figure 12. This can be stated as [29]: where K p , K i , K d are the proportional, integral, and derivative gains, m(t) is the control moment applied close to the wing root, e(t) is the error, which is tip acceleration in this case, t is the time. The upper plot of Figure 13 presents the response of the wing after store separation and focuses on wingtip oscillations for controlled and uncontrolled cases. It can be observed that the amplitude of oscillations increases without limit in the uncontrolled case, while no oscillations appear in the controlled case. The controller is active and the bottom plot of Figure 13 presents the required moment. At the point of separation, the controller takes effect and suddenly increases the amplitude of the moment to approximately 9 N m. A few steps are required to control the wing, with decreasing amplitude of control moment, which eventually vanishes.

Conclusions
This work examined the dynamic response of a high aspect ratio wing due to the separation of an externally supported store using flexible multibody dynamics. MBDyn, an open-source multibody dynamics tool, was used to model a nonlinear aeroelastic wing model, and then model the commands that release the store at a given time after the wing reaches a steady-state solution. A high aspect ratio wing, carrying a load, was used to demonstrate the benefits of using the suggested approach. The analysis showed the following: • The wing response, which is stable before the response, could be unstable after the separation of the external store. • A preliminary investigation on wing stiffness showed that a 6 percent increase in the wing structural mass is required to achieve a stable wing.

•
As an alternative to the wing stiffening, the diverging wing response after the store separation could be stabilized using a controller. In this case, a moment with maximum 9 N m amplitude is required.
This work could be followed by more complicated analysis to achieve a better understanding of the problem and to address other issues that are likely to occur after external stores are released from aircraft wings. Using the suggested method and the developed model as a basis, possible additions can be made, which are briefly mentioned, as follows: • The simulations can be repeated for a flying aircraft, including the rigid body motion, fuselage vibrational modes, and flight mechanics derivatives.

•
The store can be given aerodynamic characteristics to include the effects of store aerodynamics on the wing stability.

•
The store location, store mass, wing airfoil, and the wing structure can be optimized for a safe store separation, while ensuring that the payload carrying capacity of the aircraft is kept the same.

•
A flight envelope for stable store separation could be achieved.
Funding: Article Processing Charges were funded by Imperial College London.