Development of an efficient contact-friction model for high-fidelity cargo airdrop simulation

Leiming NING,Jichang CHEN,Mingbo TONG

College of Aerospace Engineering,Nanjing University of Aeronautics and Astronautics,Nanjing 210016,China

KEYWORDS Airdrop simulation;Cargo airdrop;Contact-friction model;Flight simulation;Multi-body dynamics;Parachute

Abstract High-fidelity cargo airdrop simulation requires the contact dynamics between an aircraft and a cargo to be modeled accurately.This paper presents a general and efficient contact-friction model for simulation of aircraft-cargo coupling dynamics during airdrops.The proposed approach has the same essence as that of the finite element node-to-segment contact formulation,which leads to a flexible,straight forward,and efficient code implementation.The formulation is developed under an arbitrary moving frame with both the aircraft and the cargo being treated as general six-degree-of-freedom rigid bodies,and thus it eliminates 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 three 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.1Since 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 a transport aircraft and airdrop systems,which offers significant economic benefits over flight tests.

The simulation of cargo airdrop has received much attention in recent years due to its highly nonlinear and strongly coupled dynamics characteristics,especially when a cargo is large.Great efforts have been devoted to modeling and simulation of cargo airdrop.For example,to study the cargo pitch motion and extraction characteristics during the aircraft extraction phase,Ke and Yang,et al.2,3modeled the cargo 2-D motion in the aircraft symmetric plane;to investigate the stability and control of an aircraft during airdrop,several authors4-6established a six-Degree-Of-Freedom(DOF)aircraft model,and treated the cargo as a 1-DOF moving particle.In these studies,Ke and Yang,et al.treated the cargo as a 3-DOF rigid body constrained in the aircraft symmetric plane and applied an extraction force computed by the parachute drag coefficient;the rest of authors either pre-defined the motion of the cargo or computed the motion by applying a constant extraction force.In all of these studies,the contact between the aircraft and the cargo was not modeled,and the unsteady coupled motion of the aircraft during the extraction phase was ignored.With regard to investigations focusing on cargo dynamics during the whole airdrop process,Cuthbert et al.7,8and Potvin et al.9reported an airdrop simulation tool named Decelerator System Simulation (DSS), which was developed by NASA and later used by US Army to support airdrop testing.DSS has the capability to simulate 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 xOz-plane of the aircraft body frame,thus only two contact points are used to detect the contact between the cargo and the aircraft,and the contact force is computed by a simple linear springdamper 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 a higher fidelity,Fraire et al.10,11established an extraction and separation model for Orion Test Vehicles in a commercial multi-body dynamics program ADAMSTM,which provides a more accurate contact loads prediction than that of its predecessor DSS.However,the aircraft motion is pre-defined by a pitch rate profile and does not response to the contact loads.More recently,Jann et al.12,13reported a PARALAB(parachute and airdrop evaluation laboratory)simulation framework developed in MATLABTM/SimulinkTMenvironment by German Aerospace Center (DLR), which has similar capabilities to those of 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 the 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 calculation of contact forces.

According to the above literature review,it can be seen that the main drawback of 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 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 simulation of coupled aircraft-cargo dynamics during airdrop.The proposed formulation treats both an aircraft and a cargo as general rigid bodies with 6 DOFs,thus the limitations of the symmetric assumption in previous studies are eliminated,and dynamics computation can be carried out in a fully-coupled manner,without the need to pre-define aircraft or cargo motion as in the semi-coupled approach mentioned above.The main contribution of our contact model and coupling algorithm is that they give FEM comparable simulation fidelity while still maintain the efficiency and flexibility for real-time simulation(which will be demonstrated by the 2nd and 3rd numerical cases in Section 6).In addition,the method aims to be flexible,efficient,and easy to implement,and thus it can be easily integrated with an existing flight simulator to add a real-time airdrop simulation capability,or it can be implemented as a standalone‘‘glue”module to bind a parachute dynamics code and a flight dynamics code together for high-fidelity multi-body cargo drop simulations.

The structure of this paper is as follows.In Section 2,we describe the contact geometries involved in cargo airdrop and the basic strategy for contact interface discretization.In Section 3,the formulations of the normal contact force and fiction force for a single contact node are firstly derived,and then the total contact force and moment under a moving contact surface frame are developed.In Section 4,the equations of motion for a general rigid body with 6 DOFs are introduced,which have been used as the state function for both an aircraft and a cargo.In addition,the algorithm for solving aircraft and cargo dynamics in a fully coupled manner is presented.Section 5 discusses algorithm implementation details and proposes some techniques for efficient code implementation.In Section 6,three numerical examples with increasing complexity are presented,of which the first two examples also serve as validation purposes.Conclusions are drawn in Section 7.

2.Contact model description

Fig.1 Typical configuration of a cargo and an aircraft roller-rail system.

Fig.2 Moving contact surface frame and contact nodes.

An essential part of contact problems is the contact geometry,and Fig.1 shows a typical configuration of a cargo and an aircraft roller-rail system.As shown in Fig.2,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 stage,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 contactfriction model has the same essence as that of the finite element node-to-segment formulation,14that is,we use one or more contact surfaces to represent the cargo,and use distributed contact nodes to represent the roller tracks and other parts fixed on the aircraft which may have contact or collision 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(for example,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.Meanwhile,the distribution of nodes can be parameterized and fully controlled by a program,without the need to manually generate a contact geometry for different aircraft configurations.

3.Contact force computation

Before computing the total contact force and moment between the aircraft and the cargo,we firstly consider a single contact node,and develop the normal contact force and tangential friction force from the contact node to the contact surface.

Instead of developing equations directly under the inertial frame,a local contact surface frame is introduced,as shown in Fig.2,which is an arbitrary moving frame attached to the contact surface,with its z-axis pointing to the inward surface normal direction.The derivation process can be greatly simplified under this local frame,and 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 fσis computed by a nonlinear spring-damper model,i.e.,

where k is the contact stiffness,c is the damping coefficient,and e is the nonlinear spring force exponent.Compared to a linear spring model with e=1.0, this nonlinear springdamper combination provides a more accurate model of the physical behavior of colliding solids.15,16As shown in Fig.2,dσis the penetration distance,i.e.,the signed distance from the contact node to the contact surface,and with the contact surface's normal vector^zcsf,it is computed by

in which rcnd/csf=rcnd-rcsfis the contact node's relative position with respect to the contact surface frame.vσ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 a velocity vcsfand an angular velocity ωcsf,vσcan be expressed as

where vcnd/csf=vcnd-vcsf-ωcsf×rcnd/csfis the contact node's relative velocity with respect to the contact surface frame.

In order to avoid a numerically non-physical‘‘pulling”force(i.e.,fσ＜0)when dσ≈0 and vσ?0,a negative contact force is trimmed to 0,namely,fσ=max(fσ,0).

It can be seen from Eq.(1)that when the contact node starts to penetrate the contact surface with a high penetration velocity(i.e.,dσ≈0 and vσ?0),fσ‘‘jumps”from zero to a high damping force discontinuously,which is both questionable physically and unfavorable numerically.To eliminate this discontinuity when dσ≈0,we can linearly increase the damping coefficient from zero to c according to the penetration distance dσ,and thus we can write c=min(cmax·dσ/δσ,cmax),where δσis the minimum penetration distance to apply the given maximum damping coefficient cmax.With the equations provided above,the final form of the contact normal force can be written in a single expression as

In Eq.(4),the contact stiffness k can be directly computed from the Young's modulus of the contact interface material(which is usually a known property for a given metallic material),while the damping coefficient c can beobtained from literature,such as Refs.14,15.

3.2.Tangential friction force

The friction force between a contact node and a contact surface is computed by a modified Coulomb friction model,which is given by

where μ is the friction coefficient.To model both the stick and sliding states in friction14,μ is computed by

where μsand μdare pre-determined 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 Eq.(5)and Eq.(6),∈fand ∈μ,which are smooth step operators for the friction force and the friction coefficient,respectively,and their values are within the range of[0,1]and determined by the contact node's motion state.

Fig.3 shows how ∈μvaries with the tangential velocity vτ.The employment of ∈μavoids convergence difficulties caused by a sudden‘‘jump”from static friction to dynamic friction,and provides a continuous and smooth numerical transition between the stick and sliding states,which corresponds to the physically observed micro slip state.17Additionally,the use of ∈μavoids switching between different algorithm branches for sticking and sliding,therefore it's more compact and efficient.

The reason for introducing ∈fin Eq.(5)is because the classical Coulomb friction model only defines the maximum static friction for the stick state,but the actual static friction force can take on any value in the interval between zero and fσμs,depending on the magnitude of the total non-frictional external force.Thus,to determine the static friction force numerically, the total non-frictional force has to be explicitly computed,which is usually difficult or even impossible in practical simulations of complex systems. Therefore, ∈fis employed to adaptively adjust the static friction force to the correct equilibrium magnitude,solely based on sliding tendency(i.e.,micro-slip velocity).

Fig.3 Smooth transition between static and dynamic frictions by ∈μ.

The smooth step operator is implemented as a third-order Hermite interpolation,which has C0and C1continuity at interval boundaries,and it can be expressed as

where u=(x-x0)/(x1-x0),and x0and x1are the left and right boundaries of the interval to be smoothly interpolated.For ∈μ=h(vτ,ξs,ξd),the velocity interval[ξs,ξd]is chosen as ξs=10-6m/s and ξd=10-3m/s,which means that when the tangent velocity vτ≤ξs,the contact node is considered as‘‘stick”state,while‘‘sliding”state when vτ≥ξd.For fτto achieve the maximum static friction force,∈fmust achieve the value of 1.0 before ∈μ,and thus ∈fis computed from∈f=h(vτ,0,κξs),where 0 ＜κ ＜1.0.

It should be noted that Eq.(7)is not the only choice for the smooth step function,and 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 transition17,18:

where vsis the Stribeck velocity,and σsis an experimentdetermined exponent.It can be seen that similar to Eq.(6),it also satisfies vτ→0 ?μ →μsand vτ?vs?μ →μd.However,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 different floating-point precisions.In addition,the present Hermite formulation has a lower computational overhead than that of 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'forces together.For simplicity,if we only consider one contact surface(i.e.,the cargo pallet bottom surface as indicated in Fig.1),then the total contact force is given by

where ncndis the number of contact nodes;andare the ith contact node's normal force and friction force,respectively;^zcsfis the contact surface frame's z-axis direction vector;is the ith contact node's tangential velocity,which is computed bySimilarly,the total contact moment vector(about the contact surface frame's origin)is given by

If more than one contact surfaces are needed for extra safety analysis,Eq.(9)and Eq.(10)need 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 the aircraft and the cargo in general motion.In this section,the focus is on the contact coupling dynamics between the aircraft and the cargo.The reference frames used in our solver and their relationships are firstly introduced.Then,equations of motion for a general 6-DOF rigid body are given,which are used as the state function of both the aircraft and the cargo.Finally,the contact coupling algorithm is detailed.

4.1.Reference frames

An important part of developing dynamic simulations is dealing with reference frames,and 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 airdrop simulation,the rotation of Earth can be safely ignored due to its negligible effects on final results.Thus,in our study,the ellipsoidal Earth is treated as an inertial system,and a non-rotating Earth-Centered,Earth-Fixed(ECEF)coordinate system E is used as the inertial frame. The World Geodetic System 1984(WGS-84)19reference ellipsoid is adopted for corresponding geodetic and gravitational computations.

Secondly,a body-fixed Front-Right-Down(FRD)coordinate system B,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 B0to denote the aircraft body frame and B1to denote the cargo body frame.

Finally,as mentioned before,for each contact surface of the cargo,a contact surface frame C is defined for contact computation.Since a 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 B1.

It should be noted that the North-East-Down(NED)coordinate system G is also used in simulation as the intermediate frame between the ECEF frame and FRD frames.If we useBXAto denote the transform matrix from frame A to frame B,the transform sequence between these frames is shown in Fig.4.Detailed definitions of ECEF,NED,and FRD frames and the transform matrix between them can be found in many modern flight dynamics textbooks such as Refs.20,21,and 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,Ev means that the velocity vector v is given in the ECEF frame.

4.2.Equations of motion

To account for the most general case,both the aircraft and the 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 approach20, the analytical mechanics approach (such as Lagrange's method and Kane's method)22,23,and the 6-D spatial vector approach.16In this study,the vector form of Newton-Euler equations is used,because it is straightforward to implement with 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:

whereBv=[u,v,w]Tis the linear velocity vector of the center of mass,andBω=[p,q,r]Tis the angular velocity vector about the center of mass;J is the 3×3 inertia tensor;BFGis the gravitational force,and since it's usually computed by a gravitational model in the ECEF frame, thus we haveBFG=BXE·EFG;BFAandBMAare the aerodynamics force and aerodynamics moment,respectively.The extra external force termBFextand moment termBMextare reserved as an interface for coupling with other components.Physically,these can be any external forces or moments from other components,for instance,for the 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 Ordinary Differential Equations(ODEs)in Eq.(11),the following kinematics ODEs are also needed:

Fig.4 Reference frames and transform sequence.

To maintain the unit norm of the rotation quaternion even in the presence of rounding errors,the method of algebraic constraint24,25is used in Eq.(12b),where kΔt ≤1(Δt is the integration time step)and λ=1-Bq·Bq.

For the convenience of code implementation,the derivatives of state variables in Eq.(11)and Eq.(12)are already arranged to the left-hand side,and thus the complete state vector for a 6-DOF body is y={Bv,Bω,Er,Bq}T.Using a vectorized representation, Eq. (11) and Eq. (12) 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 the aircraft and the cargo can be coupled together as

Eq.(15a)use the contact scheme developed in Section 3 to compute the contact force and moment,whereCFCis the total contact force acting on the contact surface,andCMCis the total contact moment on the origin of the contact surface frame.At every time step,Eq.(9)and Eq.(10)are used for the computation ofCFCandCMC,which need the motion states of contact nodes and the contact surface as inputs.The motion state of the i-th contact node(i.e.,the positionand the velocitycan be computed from the aircraft state vector as

whereEXB0is the transform matrix from the aircraft body frame to the ECEF frame,andis the constant position vector of the ith contact node with respect to the aircraft body frame. The motion state of the contact surface (i.e., the position rcsf,the velocity vcsf,and the angular velocity ωcsf)can be computed from the cargo state vector as

whereEXB1is the transform matrix from the cargo body frame to the ECEF frame,andB1rcsf/cgois the constant contact surface position vector with respect to the cargo body frame.

In Eq.(15b),B0FCandB0MCare the coupling force and moment for the aircraft,respectively.Physically,B0FCrepresents the total contact force acting on the aircraft,whileB0MCrepresents the total contact moment about the aircraft center of mass,and they are computed by

where the transform matrix from the contact surface frame to the aircraft body frame is given byB0XC=B0XE·EXB1·B1XC,and the contact surface's relative position vector with respect to the aircraft is given by

Similar to Eq.(15b),B1FCandB1MCin Eq.(15c)are the coupling force and moment for the cargo,respectively,and they are computed by

Fig.5 gives the corresponding workflow of the coupling algorithm.As shown in the figure,for every time step with a step size Δt,the aircraft ODEs,cargo ODEs,and contact model in Eq.(15)are evaluated by the ODE stepper to get the incremental state vector ΔY={Δyaft,Δycgo}T,which is then used to update the state vector Y to the next time step.

5.Numerical implementation

We implement the proposed contact-friction model as a standalone module in Python language,which can be easily integrated into existing airdrop simulation routines.Below are 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 below.

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 hand-written serial for-loops,which usually has much better performance due to the underlying utilization of parallel SIMD instructions(e.g.,SSE,AVX,etc.).Our contact formulation naturally supports vectorization,and since we have already developed the state functions in a vector form,all we need to do is using a vector array instead of a single vector in Eq.(9)and Eq.(10)for contact force computation.Programming is almost trivially easy in a language that handles array-wise and matrix operations.In our case,vectors and matrices 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 vectorized force computation,a fast range-based collision detection method is used to eliminate non-contact nodes before force computation.The collision detection method closely resembles the point-in-AABB inclusion test26,i.e.,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 ith contact node,this can be expressed in the contact surface frame as

Fig.5 Workflow of the coupling algorithm.

5.2.ODE integration

To integrate the ODEs obtained in Section 4,the odeint function from SciPy27is used,which is a wrapper around the wellknown LSODA integrator28from the FORTRAN library ODEPACK.29LSODA is the most widely distributed numerical integration method which has the capability to automatically detect ODE stiffness and switch between the non-stiff Adams method and the stiff BDF method.For every integration step,the solver automatically chooses the class of methods which is likely to be most suitable for that part of the problem,and thus,it has great reliability and efficiency regardless of the stiffness of the problem.For the detail of LSODA implementation,readers may refer to the work by Petzold.28

6.Numerical examples

To validate the above methodology and demonstrate how it can be used in cargo airdrop simulations,three numerical examples,with increasing fidelity and complexity,are presented below.

6.1.Static-dynamic friction transition

Firstly,we consider a simple yet illustrative example to demonstrate the basic usage of the above formulation and verify the contact-friction force computation by a transition between static friction of the stick state and dynamic friction of the sliding state.

As shown in Fig.6,the payload is initially sitting on a level ground,and a linearly-increasing extraction force along the global x-axis is then applied.After certain amount of time,the payload will start sliding due to the extraction force exceeding the maximum static friction force.

The input parameters used in simulation are given in Table 1.

The extraction force is given as

Fig.6 Diagram of a static-dynamic friction transition case.

Table 1 Input parameters of static-dynamic friction transition.

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 dynamic friction will be 600 N).Meanwhile,to investigate the effect of the contact node density on friction results,three simulation runs are conducted with an increasing contact node density:10 point/m,20 point/m,and 40 point/m.For a payload of 1 m in length,this leads to average 10,20,and 40 active contact nodes during simulation.Numerical integration is carried out using a fixed time step of 0.005 s.

Fig. 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 which give a 1000 N maximum static friction at t=7 s,and a 600 N average dynamic friction during sliding,which exactly match the expected analytical solutions.

During the payload sliding stage after t=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 new contact nodes ahead,resulting in a continuous contact force switching in active contact nodes.In addition,it can be observed in the zoomed inset that a higher contact node density will lead to less oscillation due to a shorter node interval.Of course,a model with a higher contact node density will be more computationally intensive,but by adopting vectorization programming and the non-contact node elimination technique given in Section 6,the present contact algorithm can be extremely efficient on modern multi-core 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 t=15 s)are 5.32 s,8.15 s,and 13.64 s,respectively.

6.2.Cargo dropped from a fixed ramp

Fig.7 Time history of the friction force.

As a second numerical example,we try to resemble a typical scenario of the cargo airdrop extraction stage,while keeping the model setup as simple as possible for validation purpose.A similar example is discussed in Ref.30.The purpose of this example is to demonstrate the integration of the proposed contact-friction model with 6-DOF rigid body dynamics,and validate simulation results by a sufficiently refined FEM model.

Fig.8 shows the basic setup of this example.A cargo is sitting on a 5 m×2 m space-fixed ramp with a 2°slope angle.Due to the static friction,the cargo is initially still.Then,a linearly-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 example problem are given in Table 2.

The moments of inertia of the cargo are given as Ixx=333.3 kg·m2and Iyy=Izz=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 commercial FEM solver LS-DYNATM,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 a 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-DYNATMbuilt-in surface-to-surface contact type31is used for the cargo-ramp contact interface,and material properties are kept the same as the friction coefficient and contact stiffness used in the present model.

Firstly,to validate the numerical implementation of 6-DOF rigid body dynamics and its integration with the proposed contact-friction model,the time histories of selected cargo state variables are given in Figs.9 and 10.

Fig.8 Diagram of a ramp drop case.

Table 2 Input parameters of cargo dropped from a fixed ramp.

Fig.9 Time history of the cargo velocity.

Fig.10 Time histories of the cargo pitch rate and pitch angle.

For translational results,Fig.9 gives the horizontal velocity(Vx)and vertical velocity(Vz)of the cargo during the whole extraction-drop process,and for rotational results,Fig.10 shows the cargo pitch rate and pitch angle during the extraction-drop process. In addition, 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 velocity and pitch motion results of the present method are in good agreement with those of the highlyrefined FEM solution during the entire process.Before the extraction(i.e.,t ＜2 s),initial oscillations can be observed in the vertical velocity and pitch rate of the FEM solution.This is because the ramp is unstrained at the start of simulation,and no initial penetration is 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 t=2 s,the influence of this oscillation is fully eliminated.

As shown in Fig.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 behavior.7Fig.11 gives the results of the normal contact force during the edge contact process.Both data curves are directly from simulation output, and no filter is 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.The high-frequency force fluctuations in the simulation results of the highly-refined FEM model are mainly caused by relative sliding and penetration between the discrete contact interfaces,and the more finite elements the contact interfaces have,the more numerical oscillation there will be.In most FEM solvers,this non-physical numerical oscillation is eliminated by applying a mild highfrequency filter to raw results.Since our approach uses one contact node per roller,thus the number of contact elements is much less than that of the FEM model,which minimizes this numerical oscillation without sacrificing much fidelity.

Fig.11 Cargo normal contact force.

6.3.Motion coupling with the aircraft

As a final example,to 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 the cargo dynamics are modeled,and the present contact-friction model is used to compute the coupling force between the aircraft and the cargo,as presented in Section 4.The components and corresponding models involved in this example are shown in Fig.12.

Fig.12 Diagram of a motion coupling case.

As shown in Fig.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 a variable added mass tensor32and shape-driven aerodynamics characteristics.1With the extraction parachute and extraction line included,an aircraft-cargo-parachute three-body coupling system with total eighteen degrees of freedom is formed.Under ideal circumstances(no crosswind,no atmospheric turbulence,etc.),this three-body eighteen-DOF model will be an overkill.However,in order to simulate many extreme cases,the fidelity of an 18-DOF model is required.For example,when under strong crosswind or turbulence,the parachute will generate a heavy and unstable side extraction load on the cargo,which will result in an unbalanced aircraft-cargo contact force and may cause safety incidents.The asymmetric contact load and motion can be correctly captured by this 18-DOF model.On the other hand,the flexibility and generality of the proposed coupling approach allows us to build a fully-coupled threebody extraction model easily,and solve the coupled dynamics efficiently.To solve the coupled dynamics of this system,the following extra equations are required for the extraction parachute:

where Eq.(23a)is the massless spring-damper line model,which calculates the tension forceEFT(given in the ECEF frame)from parachute and cargo state vectors;Eq.(23b)is the parachute state function,which is very similar to the rigid body state function given in Eq.(14),but with special mass and aerodynamics properties as mentioned above;B3denotes the parachute body frame.Since the details of the parachute dynamics are beyond the scope of this paper,for reference purpose,a typical computed extraction force profile of the 10 m2extraction parachute used in this example is directly given in Fig.13.

High-fidelity simulation and 3-D visualization are 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.

Fig.14 gives aircraft pitch motion results from the coupled simulation.It shows that the aircraft pitch rate reaches its maximum at t=5.5 s,just before the ramp clear,and the aircraft reaches the maximum 4.6°pitch angle at t=6.8 s,about 1 s after the ramp clear.

Fig.15 gives the total contact force and moment during the extraction phase.For comparison purposes,the result without a coupled aircraft motion(i.e.,the aircraft is assumed steady during the extraction phase)is also given.

Fig.13 Typical computed extraction force.

Table 3 Input parameters and simulation scenario.

Fig.14 Aircraft pitch rate and pitch angle.

Fig.15 Total contact force and total contact moment.

As shown in Fig.15,the difference between the steady and coupled results is not obvious for t ＜4.5 s,which is because the aircraft pitch angle remains unchanged before t=4.5 s,which can be observed in Fig. 14. During the interval 2.6 s ＜t ＜3.5 s,the contact moment increases rapidly due to the inflation of the extraction parachute,which is because the eccentric extraction force produces a pitch-up moment for the cargo.A sudden change of the contact moment can be observed when the cargo starts to tilt at t=5.2 s,and 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 t=5.7 s.

Two screenshots from the corresponding 3-D visualization of this simulation are given in Fig.16.For better visibility,only half of the aircraft is shown.The first screenshot is taken at t=3.7 s,when the extraction parachute is fully opened and extracting the cargo pallet through an extraction line(the extraction line is not shown);the second screenshot is taken at t=5.6 s,right before the cargo clears the edge of the ramp.

Fig.16 Screenshots of the airdrop scene visualization.

7.Conclusions

A new contact-friction modeling approach for 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 a 6-DOF cargo is described in detail.

In this paper,the method is applied to simulate three 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 is noted that results are identical to analytical solutions.The second example resembles a real cargo airdrop extraction process by extracting a cargo on a fixed ramp.This case is also modeled and simulated in a commercial FEM solver LSDYNATMfor one-on-one comparison with the present method.The critical edge contact situation before ramp clear is closely examined.The cargo motion and contact force results of the present method agree very well with those of the highlyrefined finite element model.The last example demonstrates 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 total eighteen degrees of freedom is simulated,and the capability of the present method to increase the fidelity of airdrop simulation is demonstrated.

The proposed contact modeling method should prove useful to both the flight mechanics community and the parachute dynamics community,because it provides an efficient and flexible formulation for simulation of aircraft-cargo coupling dynamics during airdrops.