Abstract
A high-fidelity cargo airdrop simulation requires the accurate modeling of the contact dynamics between an aircraft and its cargo. This paper presents a general and efficient contact-friction model for the simulation of aircraft-cargo coupling dynamics during an airdrop extraction phase. The proposed approach has the same essence as the finite element node-to-segment contact formulation, which leads to a flexible, straightforward, and efficient code implementation. The formulation is developed under an arbitrary moving frame with both aircraft and cargo treated as general six degrees-of-freedom rigid bodies, thus eliminating the restrictions of lateral symmetric assumptions in most existing methods. Moreover, the aircraft-cargo coupling algorithm is discussed in detail, and some practical implementation details are presented. The accuracy and capability of the present method are demonstrated through four numerical examples with increasing complexity and fidelity.
1. Introduction
Airdrop plays a vital role in many situations, such as military transportation, humanitarian aid, and spacecraft tests [1]. Since airdrop operations have high requirements for reliability and safety, equipment and procedures must be extensively tested and certified. Nowadays, with the rapid growth of computer processing power, airdrop simulation has become a common approach to access and evaluate the performance of the transport aircraft and airdrop systems, which offers significant economic benefits over flight tests.
The simulation of a cargo airdrop has received much attention in recent years due to its highly nonlinear and strongly coupled dynamic characteristics, especially when the cargo is large. Great efforts have been devoted to the modeling and simulation of a cargo airdrop. Although a lot of researches have been devoted to the modeling of parachute-payload multibody dynamics when the cargo is out of the aircraft [2, 3] or to studying aircraft control and stability characteristics during the cargo extraction phase [4], much less attention has been paid to accurately model the fully coupled contact dynamics between aircraft and cargo during the cargo extraction phase while the cargo is moving inside aircraft. For example, to investigate the stability and control of an aircraft during an airdrop, Feng et al. [5] modeled the aircraft as a rigid body and simplified the moving cargo as a particle with 1-D motion inside the aircraft. Xin and Shi [6] established the airdrop dynamics with crosswind and ground effect taken into account, and analyzed the closed-loop stability of the control system under weak uncertainties. Liu et al. [7, 9] also treated the moving cargo as a particle and modeled the aircraft-cargo motion by the separation body method [5]. Because the main focus of these papers is the aircraft control characteristics, to simplify the problem, in these studies the motion of the cargo is either predefined or computed by applying a constant extraction force on a particle; the contact force and moment between aircraft and cargo are not accurately modeled.
With regard to the investigations focusing on cargo dynamics, Cuthbert et al. [10, 11] and Potvin et al. [12] reported the airdrop simulation tool named DSS (decelerator system simulation), which was developed by NASA and later used by the US Army to support airdrop testing. DSS has the capability of simulating the dynamics of cargos dropped from the ramp of an aircraft, from rollout to steady decent. The contact model used in DSS assumes that cargo motion is symmetric about the xz-plane of the aircraft body frame; thus, only two contact points are used to detect the contact between the cargo and aircraft, and the contact force is computed by a simple linear spring-damper model. The motion of the aircraft in DSS is predetermined and does not result from forces and moments acting upon it. In order to model the aircraft-cargo contact dynamics with higher fidelity, Fraire et al. [13, 14] established an extraction and separation model for Orion Test Vehicles in the commercial multibody dynamics program ADAMS, which provides a more accurate contact load prediction than its predecessor DSS. However, the aircraft motion is predefined by a pitch rate profile and does not response to the contact loads. More recently, Jann et al. [15, 16] reported the PARALAB (parachute and airdrop evaluation laboratory) simulation framework developed in a MATLAB/Simulink environment by the German Aerospace Center (DLR), which has similar capabilities as DSS; besides, it also has real-time capabilities that allow the implementation and evaluation of airdrop scenarios in a cockpit simulator. The contact model of PARALAB also adopts the symmetric assumption and the two-contact-point scheme as in DSS, but a more sophisticated nonlinear spring-damper model is used for the calculation of contact forces.
According to the above literature review, it can be seen that the main drawback of the current airdrop contact modeling approaches is that they are either overly simplified that contact loads and asymmetric effects cannot be accurately predicted, or they are unduly complicated, so that a full-scale model targeting only some specific working conditions needs to be built in commercial software packages. The objective of this paper is to propose a more sophisticated and general contact modeling approach for the simulation of coupled aircraft-cargo dynamics during an airdrop. The aim of our method is to fill the gap between aircraft flight simulation and cargo airdrop simulation, with a particular focus on the high-fidelity contact coupling dynamics between aircraft and cargo. The proposed formulation treats both aircraft and cargo as general rigid bodies with 6-DOF, thus the limitations of a symmetric assumption in previous studies are eliminated, and dynamics computation can be carried out in a fully coupled manner, without the need to predefine aircraft or cargo motion like the semicoupled approach mentioned above. The main contribution of our contact model and coupling algorithm is that they give FEM-comparable simulation fidelity while still maintaining the efficiency and flexibility for real-time simulation (which will be demonstrated in the 2rd and 3rd numerical case of Section 6). In addition, the coupling algorithm aims to be compact and easy to implement; thus, it can be easily integrated with an existing flight simulator to add real-time airdrop simulation capability; or it can be implemented as a standalone “glue” module to bind the parachute dynamics code and flight dynamics code together for high-fidelity multibody cargo drop simulations.
The structure of this paper is as follows. In Section 2, we describe the contact geometries involved in a cargo airdrop and the basic strategy for contact interface discretization. Section 3 first derives the formulation of a normal contact force and fiction force for a single contact node; then, the total contact force and moment under the moving contact surface frame is developed. In Section 4, the equations of motion for a general rigid body with 6-DOF are introduced, which have been used as the state function for both aircraft and cargo. Also, the algorithm for solving aircraft and cargo dynamics in a fully coupled manner is presented. Section 5 discusses the algorithm implementation details and proposes some techniques for efficient code implementation. In Section 6, four numerical examples with increasing complexity are presented, in which the first two examples are also served for validation purposes. Conclusions are drawn in Section 7.
2. Contact Model Description
An essential part of contact problems is the contact geometry. Figure 1 shows the typical configuration of the cargo (Figure 1(a)) and aircraft roller-rail system (Figure 1(b)). As shown in Figure 1(a), the cargo consists of a packed load and a pallet, and the carrier aircraft usually has one or more restraint rails and several roller tracks. During the cargo extraction phase, the cargo will translate along the restraint rail and its pallet bottom surface will be in contact with the aircraft rollers. The basic strategy of our contact-friction model has the same essence as the finite element node-to-segment formulation [17]; that is, we use one or more contact surfaces to represent the cargo, and we use distributed contact nodes to represent the roller tracks and other parts fixed on the aircraft which may come in contact or collide with the cargo. In the simplest case, only one contact surface is needed as the pallet bottom surface, and each roller on the aircraft is treated as a contact node. If further contact safety analysis is required, the cargo can be treated as a bounding box or convex hull consisting of multiple contact surfaces, and more contact nodes can be distributed on the potential collision areas of the aircraft (e.g., the rear cargo door).
An important benefit of this node-to-segment strategy is that contact nodes can be easily adapted to different roller track configurations, and the curvature of the aircraft floor and ramp deflection (if any) can be easily modeled with properly distributed contact nodes. Also, the distribution of nodes can be parameterized and fully controlled by the program, without the need to manually generate the contact geometry for different aircraft configurations.
3. Contact Force Computation
Before computing the total contact force and moment between aircraft and cargo, we first consider a single contact node and develop the normal contact force and tangential friction force from the contact node to contact surface.
Instead of developing the equations directly under an inertial frame, a local contact surface frame is introduced, as shown in Figure 2, which is an arbitrary moving frame attached to the contact surface, with its z-axis points to the inward surface normal direction. The derivation process can be largely simplified under this local frame; the relationship between the contact surface frame and other frames will be addressed in Section 4.
3.1. Normal Contact Force
For a contact node, the normal force is computed by the following nonlinear spring-damper model: where is the contact stiffness, is the damping coefficient, and is the nonlinear spring force exponent. Compared to a linear spring model with , this nonlinear spring-damper combination provides a more accurate model of the physical behaviour of colliding solids [18, 19]. As shown in Figure 2, is the penetration distance, that is, the signed distance from the contact node to the contact surface, with the contact surface’s normal vector , which is computed by in which is the contact node’s relative position with respect to the contact surface frame. is the penetration velocity, which is the normal component of the contact node velocity relative to the contact surface. For a moving contact surface frame with velocity and angular velocity , can be expressed as where is the contact node’s relative velocity with respect to contact surface frame.
In order to avoid a numerically nonphysical “pulling” force (i.e., ) when and , the negative contact force is trimmed to 0, namely, .
It can be seen from (1) that when the contact node starts to penetrate the contact surface with a large penetration velocity (i.e., and ), “jumps” from zero to a large damping force discontinuously, which is both questionable physically and unfavorable numerically. To eliminate this discontinuity when , we can linearly increase the damping coefficient from zero to according to the penetration distance ; thus, we can write , where is the minimum penetration distance to apply the given maximum damping coefficient .
With the equations provided above, the final form of the contact normal force can be written in a single expression as
3.2. Tangential Friction Force
The friction force between a contact node and contact surface is computed by a modified Coulomb friction model, which is given by where is the friction coefficient. To model both the stick state and sliding state in friction [17], is computed by where and are predetermined static and dynamic friction coefficients, respectively.
It can be noted that, comparing with the classical Coulomb friction model, we have introduced two similar terms in (5) and (6), and , which are smooth step operators for friction force and friction coefficient, respectively. Their values are within the range of [0, 1], and determined by the contact node’s motion state.
Figure 3 shows how varies with tangential velocity . The employment of avoids the convergence difficulties caused by the sudden “jump” from static friction to dynamic friction, and provides a continuous and smooth numerical transition between the stick state and sliding state, which corresponds to the physically observed microslip state [20]. Also, the use of avoids switching between different algorithm branches for sticking and sliding; therefore, it is more compact and efficient.
The reason for introducing in (5) is because the classical Coulomb friction model only defines the maximum static friction for the stick state; the actual static friction force can take on any value in the interval between zero and , depending on the magnitude of the total nonfrictional external force. Thus, to determine the static friction force numerically, the total nonfrictional force has to be explicitly computed, which is usually difficult or even impossible in practical simulations of complex systems. Therefore, is employed to adaptively adjust the static friction force to the correct equilibrium magnitude, solely based on sliding tendency (i.e., microslip velocity).
The smooth step operator is implemented as the third order Hermite interpolation, which has and continuity at interval boundaries, and it can be expressed as where , and and are the left and right boundaries of the interval to be smoothly interpolated. For , the velocity interval is chosen as m/s and m/s, which means that when the tangent velocity is , the contact node is considered as a “stick” state, and when it is , it is considered as a “sliding” state. For to achieve the maximum static friction force, must achieve a value of 1.0 before , thus is computed from , where .
It should be noted that (7) is not the only choice for a smooth step function, it is also common to use an exponential formulation as the smooth operator in friction modeling. For instance, the Stribeck friction model uses the following expression for static-dynamic friction transition [20, 21]: where is the Stribeck velocity and is an experiment-determined exponent. It can be seen that similar to (6), it also satisfies and . But it is clear that the proposed Hermite interpolation formulation has an advantage of explicit control over the transition interval, which is very necessary for the stability of the algorithm when running under a different floating-point precision. In addition, the present Hermite formulation has lower computational overhead than the exponential formulation.
3.3. Total Contact Force and Moment
The total contact force vector from aircraft contact nodes to the cargo is obtained by summing over all contact nodes’ force together. For simplicity, if we only consider one contact surface (i.e., the cargo pallet bottom surface as indicated in Figure 1), then the total contact force is given by where is the number of contact nodes; and are the th contact node’s normal force and friction force, respectively; is the contact surface frame z-axis direction vector; and is the th contact node’s tangential velocity, and it is computed by .
Similarly, the total contact moment vector (about the contact surface frame’s origin) is given by
If more than one contact surface is needed for extra safety analysis, (9) and (10) needs to be carried out for each contact surface frame.
4. Aircraft-Cargo Coupling Dynamics
Although the derivation of the above contact formulation may seem straightforward, special care is still needed to correctly embed it into a high-fidelity flight dynamics solver with both aircraft and cargo in general motion. In this section, the focus is on the contact coupling dynamics between aircraft and cargo. The reference frames used in our solver and their relationships are first introduced. Then, equations of motion for a general 6-DOF rigid body are given, which are used as the state function of both aircraft and cargo. Finally, the contact coupling algorithm is detailed.
4.1. Reference Frames
An important part of developing dynamic simulations is dealing with reference frames. A number of frames are involved in this study.
First of all, an inertia frame is needed as the base for the derivation of dynamics equations. For an airdrop simulation, the rotation of the Earth can be safely ignored due to its negligible effects on the final results. Thus, in our study, the ellipsoidal Earth is treated as an inertial system and a nonrotating Earth-centered, Earth-fixed (ECEF) coordinate system is used as the inertial frame. The World Geodetic System 1984 (WGS-84) [22] reference ellipsoid is adopted for the corresponding geodetic and gravitational computations.
Secondly, a body-fixed front-right-down (FRD) coordinate system , with its origin at the center of body mass, is used as the body frame for 6-DOF rigid bodies. In the subsequent text, we use to denote the aircraft body frame and to denote the cargo body frame.
Finally, as mentioned before, for each contact surface of a cargo, a contact surface frame is defined for contact computation. Since the contact surface is part of the cargo, the contact surface frame can be categorized as a body frame, which has a constant transform matrix with respect to the cargo FRD frame .
It should be noted that the North-East-Down (NED) coordinate system is also used in simulation as the intermediate frame between the ECEF frame and FRD frames. If we use to denote the transform matrix from frame A to frame B, then the transform sequence between these frames appears as shown in Figure 4. The detailed definitions of the ECEF, NED, and FRD frames and the transform matrix between them can be found in many modern flight dynamics textbooks such as in [23, 24]; thus, they are not repeated here.
For the purpose of clearness, in the subsequent text we use the left superscript on a vector to specify its reference frame; for example, means that velocity vector is given in the ECEF frame.
4.2. Equations of Motion
To account for the most general case, both the aircraft and cargo are modeled as a 6-DOF rigid body. The translational and rotational dynamics of a rigid body can be formulated by multiple approaches, for instance, the Newton-Euler approach [23], analytical mechanics approach (such as Lagrange’s method and Kane’s method) [25, 26], and 6-D spatial vector approach [19]. In this study, the vector form of the Newton-Euler equations is used because it is straightforward to implement with a vectorized code and can be easily modularized in terms of program structure.
The vector form of classical Newton-Euler equations can be written as follows: where is the linear velocity vector of the center of mass, and is the angular velocity vector about the center of mass; is the 3 × 3 inertia tensor; is the gravitational force, since it is usually computed by the gravitational model in the ECEF frame; thus, we have ; and are the aerodynamics force and aerodynamics moment, respectively. The extra external force term and moment term are reserved as an interface for coupling with other components. Physically, these can be any external forces or moments from other components; for instance, for cargo this can be the tension force from the extraction line (i.e., extraction force) or the contact force from aircraft roller tracks.
To close the dynamics of the ordinary differential equations (ODEs) in (11), the following kinematics ODEs are also needed: where is the position vector of a rigid body in the ECEF frame and is the rotation quaternion of the FRD body frame with respect to the ECEF frame. The quaternion is used here instead of Euler angles because of its “all-attitude” capability and numerical advantages in a simulation [23]. is a 4 × 4 skew-symmetric matrix given by
To maintain the unit norm of the rotation quaternion even in the presence of rounding errors, the method of algebraic constraint [27, 28] is used in (14), where ( is the integration time step) and .
For the convenience of code implementation, the derivative of state variables in (11), (13), and (14) are already arranged to the left-hand side, thus the complete state vector for a 6-DOF body is .
Using a vectorized representation, (11), (13), and (14) can be generalized as where denotes the state function of a general 6-DOF rigid body.
4.3. Coupling between Moving Aircraft and Cargo
With the proposed contact model, the state functions of an aircraft and cargo can be coupled together as
Equation (17) uses the contact scheme developed in Section 3 to compute the contact force and moment, where is the total contact force acting on the contact surface and is the total contact moment on the origin of the contact surface frame. At every time step, (9) and (10) are used for the computation of and , which need the motion state of contact nodes and contact surface as input. The motion state of the th contact node (i.e., the position and velocity ) can be computed from aircraft state vector as where is the transform matrix from the aircraft body frame to the ECEF frame and is the constant position vector of the th contact node with respect to the aircraft body frame. The motion state of the contact surface (i.e., the position , velocity , and angular velocity ) can be computed from the cargo state vector as where is the transform matrix from the cargo body frame to the ECEF frame and is the constant contact surface position vector with respect to the cargo body frame.
In (18), and are the coupling force and moment for the aircraft. Physically, represents the total contact force acting on the aircraft and represents the total contact moment about the aircraft center of mass. They are computed as follows: where the transform matrix from the contact surface frame to the aircraft body frame is given by , and the contact surface relative position vector with respect to the aircraft is given by .
Similar to (18), and in (19) are the coupling force and moment for cargo. They are computed as follows:
Figure 5 gives the corresponding workflow of the coupling algorithm. As shown in the figure, for every time step with a step size , the aircraft ODEs, cargo ODEs, and contact model in (17), (18), and (19) are evaluated by the ODE stepper to get the incremental state vector , which is then used to update the state vector to the next time step.
5. Numerical Implementation
We implemented the proposed contact-friction model as a standalone module in Python language, which can be easily integrated into existing airdrop simulation routines. Section 5.1 shows some additional implementation details.
5.1. Acceleration Techniques
The main bottleneck in the proposed formulation is that the computational cost has a linear dependency on the number of contact nodes. To reduce the computational overhead with a large number of contact nodes, several acceleration techniques are used, which are described in Sections 5.1.1 and 5.1.2.
5.1.1. Vectorization
In order to compute the contact force for each node, a loop structure is needed. When using a modern high-level language (in our case, Python) for scientific computing, vectorized code is preferred over the handwritten serial for-loops, since it usually has a much better performance due to the underlying utilization of parallel SIMD instructions (e.g., SSE and AVX). Our contact formulation naturally supports vectorization. Since we have already developed the state functions in vector form, all we need to do is use a vector array instead of a single vector in (9) and (10) for contact force computation. The programming is almost trivially easy in a language that handles array-wise and matrix operations. In our case, the vectors and matrixes in equations are directly mapped into the ndarray data type from the well-known Python package NumPy.
5.1.2. Collision Detection
To further reduce the array size involved in the vectorized force computation, a fast range-based collision detection method is used to eliminate the noncontact nodes before the force computation. The collision detection method closely resembles the point-in-AABB inclusion test [29], that is, the cargo is treated as an axis-aligned bounding box (AABB), and we test whether a contact node is inside the AABB or not. For the th contact node, this can be expressed in the contact surface frame as where the boolean variable is the contact state of th contact node. Its value indicates whether the node is in contact with the cargo or not. Once the contact state of every contact node is known, the boolean array is used as the selection mask for the contact node vector array to eliminate any noncontact node. To make the algorithm even more efficient, the lazy evaluation strategy can be used in (24), (25), (26), and (27), namely, we evaluate (25) only when (24) is true, and we evaluate (27) only when both (24) and (25) are true.
5.2. ODE Integration
To integrate the ordinary differential equations (ODEs) obtained in Section 4.2, the odeint function from SciPy [30] is used, which is a wrapper around the well-known LSODA integrator [31] from the FORTRAN library ODEPACK [32]. LSODA is the most widely distributed numerical integration method which has the capability to automatically detect ODE stiffness and switch between the nonstiff Adams method and the stiff BDF method. For every integration step, the solver automatically choses the class of methods which is likely to be most suitable for that part of the problem; thus, it has great reliability and efficiency regardless of the stiffness of the problem. For the detail of the LSODA implementation, the reader may refer to the work by Petzold [31].
6. Numerical Examples
To validate the above methodology and demonstrate how it can be used in cargo airdrop simulations, four numerical examples, with increasing fidelity and complexity, are presented below.
6.1. Static-Dynamic Friction Transition
First, we consider a simple yet illustrative example to demonstrate the basic usage of the above formulation and verify the contact-friction force computation by the transition between the static friction of the stick state and the dynamic friction of the sliding state.
As shown in Figure 6, the payload is initially sitting on level ground and a linear increasing extraction force along the global -axis is then applied. After a certain amount of time, the payload will start sliding because the extraction force exceeds the maximum static friction force.
The input parameters used in the simulation are given in Table 1.
The extraction force is given as
The above parameters are intentionally chosen to be as simple as possible, so the analytical solution of this problem will be obvious (i.e., the maximum static friction will be 1000 N and the dynamic friction will be 600 N). Also, to investigate the effect of contact node density on friction results, three simulation runs are conducted with increasing contact node density: 10 pts/m, 20 pts/m, and 40 pts/m. For a payload with a length of 1 m, this leads to an average of 10, 20, and 40 active contact nodes during simulation. Numerical integration was carried out using a fixed time step of 0.005 s.
Figure 7 shows the computed total friction force during extraction. The results are essentially identical for the three cases with different contact node densities; all of them give a 1000 N maximum static friction at = 7 s and a 600 N average dynamic friction during sliding, which exactly match the expected analytical solutions.
During the payload sliding stage after = 7 s, minor force oscillations can be observed for all three cases. This is expected because during the movement, the payload will keep coming across the new contact nodes ahead, resulting in continuous contact force switching in active contact nodes. Also, it can be observed in the zooming inset that a higher contact node density will lead to less oscillation due to a shorter node interval. Of course, the model with a higher contact node density will be more computationally intensive, but by adopting vectorization programming and a noncontact node elimination technique given in Section 6, the present contact algorithm can be extremely efficient on modern multicore CPUs. To give users a rough hint, on a quad-core 3.90 GHz CPU, the average CPU time needed in three simulation runs (the numerical integration ends at = 15 s) are 5.32 s, 8.15 s, and 13.64 s, respectively.
6.2. Cargo Dropped from a Fixed Ramp
As a second numerical example, we try to resemble the typical scenario of a cargo airdrop extraction phase, while keeping the model setup as simple as possible for validation purposes. A similar example is discussed in [3]. The purpose of this example is to demonstrate the integration of the proposed contact-friction model with 6-DOF rigid body dynamics and validate the simulation results by a sufficiently refined FEM model.
Figure 8 shows the basic setup of this example. A cargo is sitting on a 5 m × 2 m space-fixed ramp with 2 slope angles. Due to the static friction, the cargo is initially still. Then, a linear increasing extraction force along the global x-axis is applied. The cargo will be extracted along the ramp until it passes the end of the ramp and starts falling down due to gravity.
The input parameters for this sample problem are given in Table 2.
The moments of the inertia of the cargo are given as = 333.3 kg·m2 and = 833.3 kg·m2. The extraction force is given as
To validate the simulation results of the present method, the problem is also modeled and solved in the commercial FEM solver LS-DYNA, which is well known for its sophisticated contact algorithms. To minimize the errors caused by a finite element mesh, a mesh refinement study is carried out and a sufficiently refined mesh has been chosen for the final comparison. In the final mesh, the cargo is modeled by 5600 solid elements, and the ramp is modeled by 3200 shell elements. The LS-DYNA built-in surface-to-surface contact type [33] is used for the cargo-ramp contact interface, and the material properties are adjusted to match the friction coefficients and contact stiffness used in the present model.
First, to validate the numerical implementation of 6-DOF rigid body dynamics and its integration with the proposed contact-friction model, the time history of selected cargo state variables are given in Figures 9 and 10.
For the translational results, Figure 9 gives the horizontal and vertical velocities of the cargo during the whole extraction-drop process. For the rotational results, Figure 10 shows the cargo pitch rate and pitch angle during the extraction-drop process. Also, to give an intuitive understanding of how the velocity and pitch motion change during the process, the state of the cargo is illustrated with indication lines and color spans in both figures. It can be seen that both the velocity and pitch motion results of the present method are in good agreement with the highly refined FEM solution during the entire process. Before the extraction (i.e., < 2 s), initial oscillations can be observed in the vertical velocity and pitch rate of the FEM solution. This is because the ramp was unstrained at the start of the simulation, and no initial penetration was imposed on the FEM model contact interface. As a result, it takes about 1 s for the FEM model to reach the balanced contact state. Since we intentionally start the extraction at < 2 s, the influence of this oscillation is fully eliminated.
As shown in Figure 8, when the cargo is partially out of the ramp, it will start to rotate due to the unbalanced pitch moment, and there will be an interval in which only the contact nodes on the ramp edge are in contact with the cargo. Accurate capturing of the contact force during this edge contact process is crucial to the prediction of the cargo rotational behaviour [10].
To validate the contact force results of the present method, Figure 11 gives the results of the normal contact force during the edge contact process. Both data curves were directly from the simulation output, and no filter was applied. Due to the highly discrete nature of the finite element method, significant force fluctuations can be observed in the FEM solution, especially during the edge contact interval. On the other hand, the present method gives a much smoother force profile due to its simplicity, and in general it shows a very good agreement with the mean value of the FEM solution.
6.3. Motion Coupling with Aircraft
To further demonstrate the practical usage of the present contact-friction model, the extraction phase of a high-fidelity cargo airdrop simulation is presented, in which both the carrier aircraft dynamics and cargo dynamics are modeled, and the present contact-friction model is used to compute the coupling force between the aircraft and cargo, as presented in Section 4. The components and corresponding models involved in this example are shown in Figure 12.
As shown in Figure 12, the extraction force is provided by the extraction line, which is connected to an extraction parachute. The tension force in the extraction line is computed by a massless spring-damper (similar to the normal contact force); the parachute is modeled as a 6-DOF rigid body with variable added mass tensor [2] and shape-driven aerodynamics characteristics [1]. With the extraction parachute and extraction line included, an aircraft-cargo-parachute three-body coupling system with a total of eighteen degrees-of-freedom is formed. To solve the coupled dynamics of this system, the following extra equations are required for the extraction parachute: where (30) is the massless spring-damper line model, which calculates the tension force (given in ECEF frame) from the parachute and cargo state vectors; (31) is the parachute state function, which is very similar to the rigid body state function given in (16), but with special mass and aerodynamic properties as mentioned above; and denotes the parachute body frame. Since the details of the parachute dynamics are beyond the scope of this paper, for reference purposes, a typical computed extraction force profile of the 10 m2 extraction parachute used in this example is directly given in Figure 13.
The high-fidelity simulation and 3-D visualization were carried out by utilizing our in-house flight simulation framework named “TRAVIS,” which was developed and maintained by the same authors and has been adopted by several industrial departments and research agencies for production use for a few years.
The key input parameters and simulation scenario are summarized in Table 3.
Since the dynamics of the 6-DOF aircraft and 6-DOF cargo are solved in a fully coupled manner, both the cargo and aircraft motion results can be computed in a single run. Figure 14 gives the aircraft pitch motion results from the coupled simulation. It shows that the aircraft pitch rate reaches its maximum at = 5.5 s, just before the ramp is clear, and the aircraft reaches the maximum 4.6 pitch angle at = 6.8 s, about 1 s after the ramp is clear.
Figure 15 gives the total contact force and moment during the extraction phase. For comparison purposes, the result without the coupled aircraft motion (i.e., the aircraft is assumed steady during the extraction phase) is also given.
As shown in Figure 15, the difference between the steady and coupled results is not obvious for < 4.5 s. This is because the aircraft pitch angle remains unchanged before = 4.5 s, which can be observed in Figure 14. During the interval at 2.6 < < 3.5 s, the contact moment increases rapidly due to the inflation of the extraction parachute. This is because the eccentric extraction force produces a pitch-up moment for the cargo. A sudden change of contact moment can be observed when the cargo starts to tilt at = 5.2 s. Then, the moment soon reaches its maximum and decreases to zero together with the contact force when the cargo clears the edge of the ramp at = 5.7 s.
Two screenshots from the corresponding 3-D visualization of this simulation are given in Figure 16. For better visibility, only half of the aircraft is shown. The first screenshot was taken at = 3.7 s; the extraction parachute is fully opened and is extracting the cargo pallet through an extraction line (the extraction line is not shown). The second screenshot was taken at = 5.6 s, right before the cargo clears the edge of the ramp.
6.4. High Extraction Load “Bounce” Effect
As a final example, in order to demonstrate the high-fidelity capability of our contact model, we try to reproduce the “bounce” effect in real airdrops when the cargo is extracted by a large parachute force. It is not possible for the traditional contact treatments to capture this phenomenon in a fully 3-D manner because the cargo is either simplified as a particle or assumed to be moving along the 2-D aircraft floor in these models. The setup of this simulation is basically the same as the last case, except that the cargo mass is reduced to 5000 kg and the nominal area of the extraction parachute is increased to 14 m2, and the aircraft ramp deflection angle is 5. Other inputs are kept the same as Table 3.
Figure 17 gives a screenshot of the real-time simulation of this case, in which the cargo trajectory is drawn. The “bounce” effect can be clearly observed in the figure, which is caused by the off-track motion of the cargo when a large extraction force is applied. The corresponding contact force and moment is given in Figure 18.
It can be seen in Figure 18 that a large recontact load occurs when the cargo bounces on the ramp, which may bring a major security risk to both aircraft and cargo. Thus, this “bounce” effect needs to be eliminated in the design stage or later optimization stage, and our contact model provides sufficient fidelity and efficiency to capture this phenomenon in a real-time simulation.
7. Conclusions
A new contact-friction modeling approach for the simulation of aircraft-cargo coupling dynamics during an airdrop is presented in this paper, and the coupling scheme between a general 6-DOF aircraft and 6-DOF cargo is described in detail.
In this paper, the method is applied to simulate four numerical examples with increasing complexity. In the first example, a simple static-dynamic friction transition process is simulated to validate the friction force computation scheme, and it was noted that the results are identical to the analytical solutions. The second example resembles the real cargo airdrop extraction process by extracting the cargo on a fixed ramp. This case was also modeled and simulated in the commercial FEM solver LS-DYNA for one-on-one comparison with the present method. The critical edge contact situation before the ramp was clear was closely examined. The cargo motion and contact force results of the present method agrees very well with the highly refined finite element model. The last example demonstrated the practical usage of the present contact-friction model in a high-fidelity cargo airdrop simulation. In this case, the dynamics of an aircraft-cargo-parachute three-body coupling system with a total of eighteen degrees-of-freedom was simulated, and the capability of the present method to increase the fidelity of the airdrop simulation was demonstrated.
The proposed contact modeling method should prove useful to both the flight mechanics community and parachute dynamics community, because it provides an efficient and flexible formulation for the simulation of aircraft-cargo coupling dynamics during airdrops.
Data Availability
The simulation result data used to support the findings of this study are included within the files of the supplementary materials.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Supplementary Materials
The simulation data files of the above numerical examples are provided in the supplementary materials, along with all of the full resolution figure files in this paper. (Supplementary Materials)