Multibody dynamics analysis of the human upper body for rotorcraft–pilot interaction

The study of the biodynamic response of helicopter passengers and pilots, when excited by rotorcraft vibrations that are transmitted through the seat and, for the latter, the control inceptors, is of great importance in different areas of aircraft design. Handling qualities are affected by the proneness of the aircraft to give rise to adverse interactions, an unwanted quality that can be captured by the so-called biodynamic feedthrough. On the other hand, the transmissibility of vibrations, especially from the seat to the head, affects the comfort of pilots and passengers during flight. Detailed and parametrised multibody modelling of the human upper body can provide a strong base to support design decisions justified by a first-principles approach. In this work, a multibody model of the upper body is formed by connecting a previously developed detailed model of the arms to a similarly detailed model of the spine. The whole model can be adapted to a specific subject, identified by age, gender, weight and height. The spine model and the scaling procedure have been validated using the experimental results for seat to head transmissibility. The coupled spine-arms model is used to evaluate the biodynamic response in terms of involuntary motion induced on the control inceptors, including the related nonlinearities.


Introduction
Rotorcraft represents an intrinsically harsh vibration environment, due to their inherent structural and dynamical characteristics. Pilots' and passengers' bodies are, therefore, generally subject to a significant vibratory excitation during flight, with twofold consequences: on the one side, both pilots and passengers can feel discomfort and, in severe cases, also suffer healthrelated issues [18][19][20]; on the other hand, in the case of pilots, discomfort is also associated with a degradation in the efficiency in performing a task. Furthermore, their upper body biomechanical response is transmitted to the command inceptors, a phenomenon that goes under the name of biodynamic feedthrough [35,36,38].
Both discomfort and biodynamic feedthrough depend on the human body vibratory response: the efficient estimation of the latter is therefore of great importance, especially when in early aircraft design stage. This aspect is particularly evident when the wide variability of the parameters that influence the biomechanical behaviour of the human body are taken into account: mechanical properties and geometry above all that further depend on anthropometric parameters like age, gender, weight and stature. The confidence with which design choices are made in the early design stages is directly related to the robustness with which estimation of the vehicle and its occupants dynamical behaviour can be carried out.
First-principles approaches, like multibody modelling, are to be preferred, in this context, precisely for their ability to account for all the different sources of variability in physics-based fashion. Fully parametrised biomechanical multibody models of the human body can support the evaluation of the fitness of design choices with respect to comfort, and the robustness of the design with respect to the insurgence of adverse interaction phenomena.
This work presents the current status of an enduring research effort by the authors' group, specifically focused on the multibody modelling of the human upper body for the purpose of enhancing the comprehension of the human-machine interaction between the rotorcraft and its occupants. A complete model able to comprehend the coupled behaviour of the spine and the upper limbs has been developed using the free generalpurpose multibody solver MBDyn 1 [27], capable of direct and inverse dynamics analysis, and of reduced order model (ROM) extraction [40]. The spine model is the main focus of the presented work: it has been developed incrementally in the last several years [32] and received major focus towards its thorough validation, and its conjunction with the upper limbs model, also developed using MBDyn that was originally presented in [30].

Human spine biomechanical models
Approaches to the numerical modelling of the human body dynamical behaviour, in the context of evaluating human-machine interactions in both comfort and handling qualities-related aspects vary significantly in previous experiences reported in the literature. They can, nonetheless, be broadly divided into three categories: Lumped parameters models (LPM), finite element models (FEM) and multibody dynamics (MBD) models. 1 http://www.mbdyn.org/ 2.1 Lumped parameter models Lumped parameter models make use of elementary mechanical elements, such as lumped masses and linear viscoelastic elements. The main advantage of these models is their low computational cost and ease of parameter identification. The core of LPM is the identification of the human body as a simplified mechanical system. The parameters are tuned to fit the biomechanical characteristics of a specific group of subjects. On the other hand, one disadvantage of LPM is that they are not consistent with the actual parameters of human anatomy and biodynamics, so they may not fully reflect the response of each part of the human body and are intrinsically difficult, or impossible, to extend outside the domain of parameters values upon which they are developed and validated. Moreover, they are often used to describe a uni-directional dynamic response of the human body because of compact expression and effective performance of the models only requiring as parameters the equivalent mass, stiffness and damping.
One of the first lumped models in this context was proposed in [10]. Starting from the measurement of mechanical impedance of eight people with different heights, weights and ages, a one degree-of-freedom (DOF) linear model was proposed. In [43], a one-DOF nonlinear model was proposed and applied to the analysis of the dynamic response of human bodies during helicopter landing. Reference [51] presented a four-DOF series-to-parallel linear model (shown in Fig. 1c) to describe the dynamic response of seated occupants subjected to vibrations. This nonlinear model was optimised in [1] using a genetic algorithm (GA). In [6], a measure of the vertical driving-point mechanical impedance of seated vehicle drivers in the range between 0.625-10 Hz was proposed; a four-DOF linear model was derived, as shown in Fig. 1d, whose parameters were estimated to satisfy both the measured driving-point mechanical impedance and the seat-tohead transmissibility (STHT). The four degrees of freedom are associated with the vertical displacement of the masses of head and neck, upper torso, viscerae and lower torso. An optimisation of the model presented in [6] was done in [56] using an average-weighted genetic algorithm. Several other lumped parameters models were proposed: the interested reader is referred to Ref. [23] for a thorough review of the subject.
In [3], a methodology is proposed for systematically identifying the best configuration or structure of Fig. 1 Lumped parameters models: a, b are the 4DO F 12−6 and 4DO F 14−9 of Ref. [3]; c is the model of Refs. [1,51]; d is that of Refs. [6,56]. The model elements can be generally associated with the structures of head and neck, upper torso, viscerae and lower torso; however, in some cases they do not have a strict physical meaning and therefore cannot be associated with specific anatomical structures a four-DOF human vibration model and for its parameter identification. The results of this optimisation are shown in Fig. 1a, b. The models were calibrated using the frequency response functions recommended by ISO 5982 [15].

Finite element models
Finite element models are more complete, with respect to LPM ones. First of all, they respect the geometry of the human body; furthermore, they can be parametrised, to represent individuals with different properties. Moreover, the output of the analysis is far richer: for example, in a frequency response analysis the acceleration and the internal load of each vertebra of the spine can be measured. The two major drawbacks are the higher computational cost and the difficult identification of the mechanical parameters of the model, compared to LPMs. FEM models can be divided into two categories: discrete and continuous. The discrete models treat the spine as a structure made by rigid elements, representing the vertebral bodies that are connected through deformable elements representing the intervertebral discs. On the contrary, the continuum models treat the spine as a homogeneous beam.
In [47], each vertebra was modelled separately in an eight-DOF nonlinear discrete model. However, the model is restricted to only the axial direction. In [34], a discrete parameter model of the human body under a variety of impact situations was developed; in detail, the pilot ejection problem was deeply investigated. These two models were restricted to one-or two-dimensional behaviour and did not consider the interaction of the spine with other parts of the torso, such as the rib cage and the viscerae.
In [5], a three-dimensional, discrete mathematical model of the human spine, torso and head was developed. This model is based on a small strain, large displacement formulation. In [4], the visceral masses, dampings and stiffnesses and the visceral modelling, were tuned in order to better match the experimental results. The visceral wall can, in fact, be viewed as a secondary path of force transmission between the seat and the upper torso: in the proposed model, forces are transmitted through the rigid link that connects the topmost visceral element to T10. The model presented in [21], an evolution of the one proposed in [4], is a 2D FEM model, which allows movements only in the sagittal plane. It models the spine, viscera, head, pelvis and buttock tissues, using beam, spring and mass elements. The geometrical and mechanical parameters are based on those presented in [4], but some of them have been modified to match the mode shapes obtained by the experimental work [22]. The model is entirely linear; it includes 134 elements and 87 master degrees of freedom.
In [46], a spine FEM model inspired by that of [21] was proposed by the authors (Cf. Fig. 2), for the purpose of rotorcraft comfort assessment: with respect to  [46] the model presented in [21], the proposed model can capture the full 3D behaviour of the spine about a nominal reference condition.

Multibody models
Multibody dynamics (MBD) adds flexibility to LPM with the ease of constraint formulation. Furthermore, flexible elements include and extend the capabilities of FEM models. The great advantage of MBD is that it can capture effects related to nonlinearities, especially the ones originating from 3D geometry.
In [53], a multibody model of the spine having ten DOF was presented: this model was restricted to the sagittal plane; the cervical and thoracic sections were expressed as one rigid body.
In [12], a two-dimensional 20-DOF model that can be used to study the vertical and fore-aft vibrations was proposed. The parameters were optimised using a GA algorithm.
A full multibody model of the spine was proposed in [49] for an accurate assessment of seated body vibration, but also this model was limited in the sagittal plane. In [50], the model was improved in order to be able to capture the three-dimensional motion of the spine.
The advantages and disadvantages of each type of models are summarised in Table 1.

Human upper body model
The proposed multibody model of the upper body consists of a detailed model of the torso, of the upper limbs, and their interface with the surrounding environment, typically an aircraft or rotorcraft cockpit.

Torso model
The torso model comprises 34 rigid bodies associated with nodes (i.e. entities possessing degrees of freedom) placed in correspondence of vertebrae from S1 to C1, the head, and of 8 visceral masses elastically connected to vertebrae from S1 to T10. The model follows the concept originally proposed in 1997 by Kitazaki and Griffin [21], in turn based on the database provided in 1978 by Privitzer and Belytschko [4].
Each vertebral section node is located in the centre of the corresponding vertebra's body, with axis z aligned with the local tangent to the curve described by the spine's longitudinal axis, axis y directed laterally and axis x pointing anteriorly (Fig. 4).
The spine itself is composed of 25 vertebral nodes connected by 24 viscoelastic 3D deformable elements, acting on both the relative displacements and rotations between adjacent nodes. They are represented in a simplified way in Fig. 3 by the linear and angular springs of elastic constants k a and k θ . Linear viscoelastic constitutive laws relate the strains of the deformable elements to their internal forces and moments. The head and the Sacrum (S1) are modelled as rigid bodies and are, respectively, connected to vertebrae C1 and L5 by linear viscoelastic elements.
When the analysis is restricted to the sagittal plane, vertebral nodes are connected to each other by algebraic constraints limiting their relative degrees of freedom to sliding along the spine axis and rotating about the lateral axis. Otherwise, when the model is used for three-dimensional simulations, all of the relative rotation degrees of freedom of the vertebral nodes are unconstrained, whereas the translation constraints remain in place (Cf. Fig. 3). The intervertebral constraints mimic the behaviour of the anatomical joints: the superior and inferior faces of the vertebrae bodies are joined by a symphysis joint, allowing all the relative rotations and relative displacements permitted by  the deformations of the fibrocartilaginous disc that lies between them; however, the intervertebral relative displacement in the transverse plane is constrained by a Fig. 4 Each node related to a vertebral section is located in the centre of the corresponding vertebra's body, with axis z i along the local tangent to the spine axis curve and axis y i pointing laterally facet joint, a synovial joint that acts on the posterior processes.
The upper and lower visceral masses are rigidly connected to vertebrae T10 and S1, while the other ones are connected to the respective vertebra node, in the section between T11 and L5, by kinematic constraints allowing only their relative displacement in the sagittal plane. Linear viscoelastic elements act on the relative viscera-vertebra displacements along the local x i and z i directions (Cf. Figs. 3, 4). Visceral nodes do not possess rotational degrees of freedom: as such, the visceral masses are considered as point masses. The elastic elements are depicted as linear springs in Fig. 3: k vs denotes the elastic constants of the springs connecting a vertebral node and a visceral node, while k ss denotes the elastic constants of the springs connecting two visceral nodes.
The detailed modelling of the muscular structures of the torso has been avoided, since the additional modelling effort would have not been justified for the application cases for which the model is currently utilised for. Torso muscles of both rotorcraft pilots and passengers are involved in postural control, i.e. typically they act to maintain a steady reference condition. The fully detailed modelling would allow to greatly extend the range of validity of the model, especially with regard to poses deviating considerably from reference one (sitting rest position). It would, however, also make the model underdetermined due to overactuation, and special techniques would have to be applied to estimate the equivalent impedance at each joint. While similar experiences have been successfully carried out by the authors, regarding, in particular, the upper limbs' impedance at the control inceptors [28,29], it was concluded that for this particular case the simpler approach relying on lumped viscoelastic elements is sufficient for the targeted analyses.
The pelvic region is modelled taking into account the compliance of the buttocks' tissue, using two viscoelastic elements that connect the Sacrum with a node representing the mean interface point between the buttocks and the resting surface. A single rigid body, accounting also for a third of the mass of the tights, in addition to that of the buttocks and pelvic structures, is connected to the pelvis node.
Geometry and inertial parameters are adapted to represent a generic subject possessing the desired anthropometric characteristics of age, gender, stature and weight, as detailed in the next section. Table 4 lists the stiffness and damping values for a reference subject of 34 years, 84 kg and 1.783 m. The values K a and C a are associated with the axial springs and dampers that act in z direction; K b i and C b i are associated with rotational springs and dampers that act about axis i. Table 5 lists the mass value, the moments of inertia and the position of each node for the reference subject.
The total number of degrees of freedom of the system, before constraints are enforced, is 228. After constraint enforcement, 103 degrees of freedom remain.

Upper limbs model
The multibody spine model of the spine has been joined with the upper limbs multibody model developed in [30,31] based on [37].
The model comprises four nodes (and the associated rigid bodies) per limb, associated with the humerus, radius, ulna and hand. The total number of unconstrained degrees of freedom is thus 24. The hand is condensed into a single rigid body since the applications that are simulated usually involve gripping tasks, in which the dynamics of the fingers can be disregarded. The nodes are connected by ideal algebraic constraints that model the following articular joints: The application of the algebraic constraints eliminates 17 degrees of freedom, resulting in a 7 degrees of freedom system. Each limb is actuated by 25 muscle fascicles, modelled by active viscoelastic rod elements. The force F i exerted by the ith muscle therefore depends on the distance between the origin and insertion points, rigidly attached to the respective nodes, on its deriva-tive˙ and on the non-dimensional activation coefficient a, ranging from 0 to 1: where the definitions of f 1 ( ), f 2 (˙ ) and f 3 ( ) follow the simplified Hill-type model [52] presented in [37].
Since the applications shown in this work involve tasks in which the elevation angle of the humerus is consistently kept at very low values, the shoulder complex is not modelled in detail, disregarding the clavicle and the scapula. This is the case of piloting tasks: the expected effect of the scapula and clavicle motion on the shoulder kinematics is very limited. The glenohumeral joint is represented, in the simplified scheme, by a spherical joint located at the glenoid fossa. The interested reader is referred to the cited literature for further details regarding the upper limbs model. The complete upper body model has been coupled with the seat model proposed in [9]. The seat is modelled as a single rigid body, suspended to the airframe floor through a linear viscoelastic element, tuned in order to limit the transfer of loads to the helicopter seat occupant. The cushion is modelled as a unilateral, nonlinear viscoelastic element having the following constitutive law: where a and b are positive real coefficients and is the deformation of the viscoelastic element. This unilateral constitutive law has been chosen in order to have null viscoelastic force when the element is undergoing positive (tensile) straining and an exponentially increasing reaction force to negative (compressive) deformation. The backrest has been modelled using five unilateral viscoelastic elements, which link the vertebrae from T4 to T8 to the backrest node. These elements share the same constitutive law with the cushion. The backrest node is rigidly connected to the seat node. A schematic representation of the seat model is shown in Fig. 5.

Model parametrisation
The model is mainly intended for the estimation of the vibratory response of the upper body of rotorcraft passengers and pilots, thus focusing on dynamics in the 1-30 Hz frequency range. In this context, the model should be able to reproduce the average, or most plausible subject possessing a certain set of anthropometric characteristics, while it is not intended to be purely subject-specific. The parametrisation is based on an anthropometric dataset s that comprises age a, body mass index BMI, stature h and gender g: The statistical model proposed by Shi et al. [42] is used to generate the most plausible rib cage geometry associated with the anthropometric parameters. The model was established analysing the thorax computed tomography scans of 89 subjects, divided into eight age groups and of both sexes. Thanks to image segmentation the position of 464 landmarks were identified on the right side of each subject's rib cage (Fig. 6).
The correlations between the landmarks positions and the anthropometric parameters were studied through a principal component analysis (PCA), following the procedure outlined by Allen et al. [2]. After generating the rib cage landmarks, a bounding box is fitted to them, and its dimensions x, y, z are compared to those of the reference subject of the work of Kitazaki and Griffin [21], x 0 , y 0 , z 0 , to yield three scaling coefficients: Since the anthropometric characteristics of the modelled subject are not reported in [21], they were estimated a posteriori using a sequential quadratic programming optimisation algorithm find the set s 0 that minimises the mean squared error of the distances between the estimated landmarks that lie closer to the vertebrae transverse processes and their corresponding locations estimated from the data of Kitazaki and Griffin. The resulting reference subject is a 34-yearold male, 1.78 m tall and weighting 84 kg, thus having a BMI of approximately 26.5. Using as reference the normal configuration found in [21], the initial positions of the nodes for the mod-

Mass scaling
Each rigid body of the upper limbs model, comprising the bone and the associated soft tissue, is scaled with the total weight using linear regression relationships identified in large survey studies [7,8] (Cf. Appendix 1). According to [11], the mass of the trunk M t of the subject accounts for the 68% of the total body mass M b , excluding the upper limbs mass M l : The masses of the viscerae M v are estimated as 20% of the total body mass, excluding the upper limbs: Then, the ratio between each trunk and visceral mass element and the total trunk and viscerae mass is defined as: where the values provided in [21] are used for M t i and M v i . The scaled value of each mass element is then obtained by multiplying the mass ratios T i and V i by the total trunk and viscerae mass.

Scaling of stiffness and damping coefficients
The parameters of the intervertebral elements' constitutive laws are scaled as follows: 1. the ratio between each stiffness or damping coefficient and the corresponding mass element is evaluated, using the parameters of the model in [21]; 2. the value of this ratio is then multiplied by the corresponding mass scaled for the specific subject; 3. the resulting stiffness coefficients are modified using a second-order polynomial P(BMI) (Cf. Appendix 1) that correlates the BMI with the first resonance frequency, and the damping factor that corresponds to this peak, where k i KG and m i KG are the stiffness and mass parameters of the Kitazaki and Griffin model [21].

Solution strategies
The multibody models require specialised solution strategies to cope with their different levels of indeterminacy. Indeed, the upper limbs model is both kinematically underdetermined and overactuated. The spine model is also kinematically underdetermined [30,55]. The kinematic underdeterminacy of the upper limbs model is evidenced by the fact that-referring to a single limb-when the motion of the hand is full prescribed, one degree of freedom still remains: namely, it is possible to perform the same task with different angles of humerus elevation. The spine model presents an even higher degree of underdeterminacy if only the motion of the head is prescribed: the same head trajectory is obtainable by a wide variety of spine con-figurations, thanks to the comparably high number of allowed relative vertebral configurations. The upper limb model overactuation is, on the contrary, due to fact that usually multiple muscles act about the same articular joint, in numerous agonist/antagonist pairs. The total torque produced about the articular joint is therefore a result of the central nervous system (CNS) actuation strategy that usually leverages muscular synergies. It is usually postulated [13,54] that the central nervous (CNS) acts in order to optimise the energy expenditure towards the achievement of different targets: minimisation of the total activation, minimisation of the metabolic cost and maximisation of the limb equivalent impedance are some possible examples.
For the typical applications, i.e. the estimation of the vibratory response of helicopter pilots and occupants for comfort assessment and the evaluation of proneness to rotorcraft-pilot couplings of existing and concept aircraft layouts, the general solution strategy is therefore made of several, cascaded, simulation steps: 1. from the desired anthropometric dataset, the model geometry and structural parameters are generated according to the procedure outlined in Sect. 4; 2. an underdetermined inverse kinematics analysis is performed, in which both the pose of the spine and of the limbs are estimated imposing the motion to the hands and to the head; 3. an inverse dynamics analysis problem is solved, to yield estimates of the activations of the upper limbs muscles; 4. a direct dynamics analysis, aimed at estimating the effects of the active and reflexive muscular forces in the upper limbs, and muscular intervertebral moments, is carried out; 5. as an optional step, a generalised eigenanalysis, directly performed on differential-algebraic equations (DAE) system, in descriptor form, to extract reduced order models (ROMs).
The steps 2-5 will be briefly described hereinafter, highlighting the most important aspects and referring to the related literature for further information.

Inverse kinematics
To obtain a square problem in the kinematic inversion when computing the pose of the system, a series of static problems is set up, in which ergonomy springs act on the redundant degrees of freedom [16]. The stiffness coefficients of the springs act as penalty coefficients for the motion of the degrees of freedom they are connected to. For example, they can be crafted to minimise the norm of the internal bending moment in the sagittal plane due to the weight, or the norm of the deviation with respect to a nominal reference condition. The inverse kinematics problem can be stated, following the approach outlined in [17], directly at the configuration level, as the problem of finding the static equilibrium configuration of the system, augmented by a set of "dummy", or ergonomy, springs, respecting the imposed motion and the algebraic constraints. Therefore, the equivalent problem can be stated as the constrained minimisation of the springs' elastic energy yielding where the algebraic constraints equations are φ(x) = 0, the imposed displacements ψ(x) = α(t) and θ (x) represents the articular joint angles. The stiffness K(θ ) can be crafted to penalise configurations in which the angles differ significantly from the intermediate ones: the choice is made on the assumption that the intermediate configurations are the most ergonomic [30]. A high degree even polynomial function of θ , centred around θ ergo is used. In the case of the spine model, the same intervertebral viscoelastic elements that are described in Sect. 3 are also used as ergonomy elements. Figure 7 shows three successive configuration of the vertebral nodes during an inverse kinematics analysis. In particular, the initial, intermediate and final configuration are shown: initially, the spine is assembled using the reference coordinates of Table 1. The head is then rigidly rotated in order to align the local z-axis of the associated node with respect to the global z-axis. When the desired head configuration is reached, the straining of the intervertebral elements is stored. In subsequent analyses, it is added to the same elements as the value of prestrain.
The solution is completed at the velocity and acceleration levels by solving other two constrained optimisation problems, in which the weighted distance of the solution from the reference one x ref , obtained by numerical differentiation of the one at configuration level, is minimised. In both cases, the mass matrix M is used for weighting. At the velocity level, the problem is therefore stated as the minimisation of with respect toẋ, λ , and μ , yielding The acceleration problem is analogous.

Upper limbs inverse dynamics and muscular activations
Once the inverse kinematics problem has been solved, the total resulting muscular torques about the articular joints of the upper limbs can be computed through an inverse dynamic analysis. This step is needed to cope with the upper limbs overactuation: each joint torque is produced by multiple muscles, usually acting together. The problem is fully determined, in this case, as it can be stated on finding the torques c and the reaction forces' Lagrange multipliers λ that, at each instant of time of interest, are able to produce the desired motion x(t),ẋ(t),ẍ(t): yielding the linear and square problem (15) where θ /x is the Jacobian of the joint coordinates θ(x) with respect to the generalised coordinates x.
When the upper limbs joint torques are known, the muscular activations required to produce them can be computed. The problem is again underdetermined: where B is the matrix of the muscle forces' moment arms with respect to the nodes, and (·) + denotes pseudo-inversion. The muscle forces, as clearly shown by Eq. (1), are composed by an active and a passive part. The diagonal matrix F 0 contains the peak isometric contraction forces of the upper limb muscles. Matrix F 12 , also diagonal, contains the products of functions f 1 ( i ) and f 2 (˙ i ). They, respectively, represent the ratio, with respect to its peak isometric force, between the maximum force that the ith muscle can produce at the current length i and elongation veloc-ity˙ i . Matrix F 3 is diagonal as well; its elements are the values of the f 3 ( i ) function, representing the ratio between the passive force of the ith muscle, including the contribution of the tendon, at the current length i , and the peak isometric contraction force. A formal solution of the problem could be found by pseudo-inverting matrix (θ + /x ) T BF 0 F 12 ; however, such solution would not guarantee the resulting activations to be limited in the admissible range 0 ≤ a i ≤ 1, necessary to ensure that the force produced by each muscle is purely contractile and does not exceed the maximum value at the current length and contraction velocity, as dictated by the nature of the skeletal muscles mechanics [52]. Therefore, an optimisation problem formulated as the minimisation of a cost function J (a) that can be adapted to fit the specific simulated scenario is instead solved, subjected to the constraints of satisfying the equivalence of Eq. (16), and the activation bounds.
Another consequence of the indeterminacy of the problem is that the computed activation level can be augmented by any contribution that does not alter the total joint torques. Such contributions have been labelled by the authors torque-less activation modes (TLAMs): they can be useful, for example, to modulate the limbs equivalent impedance at the control incep-  [30], as a function of the pilot estimated workload or stress state.

Direct dynamics
When the activations in the upper limb muscles are known, a direct analysis of the desired scenario can be run. The intervertebral stiffness and damping coefficients of Table 4 already take into account the active (i.e. caused by proprioception) muscular impedance and are kept constant throughout the simulations. The latter assumption can be justified noting that the deviation from the reference condition introduced by the prestrain estimated in the inverse kinematics step is generally limited for the applications the model has been designed for (Cf. Fig. 7).
The upper limb muscles activations estimated in the inverse dynamics step, and the subsequent optimisation problem, are augmented during the direct dynamics phase with a contribution proportional the muscles' length variation and elongation velocity, with respect to the reference values estimated during the inverse kinematics step. For the single muscle this contribution, which models the reflexive part of the muscular activation behaviour [28,45], is (17) where ref is the reference muscle length in the pose estimated in the inverse kinematics step and k p ,k d are proportional and derivative gains. The total activation of each muscle bundle in the upper limb therefore is where a 0 is the reference activation, a r the reflexive part and a linear combination b of TLAMs is introduced, scaled by the gain coefficients K TLAM .

ROM extraction
The extraction of reduced order models is based on the direct linearisation of the full DAE problem [26] Mẋ = ṗ where the first block row is the definition of the momenta p, the second block row represents the generalised equilibria, with f = f(ẋ, x, t) + (Mẋ) /x and f is applied forces vector. The third block row contains the algebraic constraints' equations. The linearisation yields the descriptor form state space system with K = Mẍ − f /x and D = −f /ẋ Eigenvalues and eigenvectors of the regular matrix pencil characterising the homogeneous solution of the (20) can be found, for example, through a generalised Schur, or QZ, decomposition [33,44] of the pencil matrices. The system can then be reduced to ODE form, eliminating the algebraic state variables in a staggered manner. More details can be found in [41].

Comparison with experimental results
The biodynamic response of seated human subjects to whole body vibration of different type and magnitude has been widely investigated, frequently together In [14], the measurements results regarding the vertical apparent mass (AM) of 60 seated subject, including 24 men, 24 women and 12 children, are presented. The measurements were taken using as input a random acceleration at frequencies up to 20 Hz. In [6], the driving-point mechanical impedance of 7 male subjects was measured. The subjects were seated on a rigid seat and subjected to ten different acceleration excitations. In [22], eight subjects were exposed to vertical random vibration while maintaining three different postures on a rigid seat without backrest. The modes in the frequency bandwidth under 10 Hz were extracted. The principal resonance of the human upper body was observed at about 5 Hz. These experimental results were used to validate the FEM model presented in [21].
In [48], the AM of 80 seated adults (41 males and 39 females) were measured at frequencies between  [24] and numerical results from Valentini [50] 0.6-20 Hz. In [25], the AM and STHT were investigated by exposing the body to vibrations applied in the fore-aft, lateral and vertical direction. The ISO [15] defined the ideal ranges of AM and STHT magnitude and phase responses for subject seated without back support and exposed to vertical vibration of magnitude up to 5 m s −2 . A synthesis of the reported data regarding the biodynamic responses of the human body exposed to vibrations along different directions and the associated experimental conditions can be found in [39].
The validity of the scaling procedure of the presented model has been tested by comparing the experimental apparent mass between 0-20 Hz of three subjects reported in [48].
The corresponding anthropometric parameters are reported in Table 2.
The apparent mass is calculated from the ratio of the cross-spectral density between the force F zs and acceleration at the seat A z s , to the power spectral density of the acceleration at the seat. The H 1 estimator is obtained: where the exponent (·) * indicates complex conjugation. The three experimental and numerical apparent mass curves are compared in Fig. 8. The model correlates well with the experimental results. Both the static mass and the first resonance frequency are well captured.
What is here proposed is intended as a full threedimensional nonlinear model of the upper body; as a consequence, to validate its behaviour along the three axes, the seat-to-head transmissibility (STHT) up to a frequency of 20 Hz has been computed and compared with experimental data from [25] and with the results from the multibody model proposed in [50]. The STHT is computed as the ratio between the cross-spectral density of the accelerations at the head A h i and at the buttock A b j and the power-spectral density of the acceleration at the buttock: The results are shown in Fig. 9. The model correlates relatively well with the experimental results for frequencies up to 10 Hz with respect to same-axis excita-tion and measurement, whereas non-negligible differences can be observed in the off-axis term H xz , i.e. in the head longitudinal response to a vertical acceleration input. Work is currently being done by the authors to understand the origin of the discrepancy that is probably due to the distribution of stiffness and damping coefficients in the intervertebral elements' constitutive laws. It should be noted, however, that the off-axis response of the human upper body has been shown to be characterised by a relatively high level of variability [57].

Biodynamic feedthrough
The complete upper body model of the reference subject has been used to perform a frequency sweep analysis in the frequency range of 0.5-12.5 Hz, evaluating the frequency response function (FRF) between the cockpit floor vertical acceleration and the rotation of the collective and cyclic levers, usually termed biodynamic feedthrough (BDFT). In Eq. (23), θ i is the ith control input (collective lever rotation, cyclic lever rotation in longitudinal or lateral direction), θ 0 is the trim angle of the lever, andz g is the imposed floor acceleration in the vertical direction: The plots in Fig. 10 compare the FRFs of the complete upper body model with those of a model consisting of only the upper limbs, for a vertical acceleration input of 0.2 g. Table 3 presents the values of frequency and amplitude at resonance for the two models. As one can see from the three FRFs, when the complete dynamics of the spine is considered the magnitude of the resonance peak increases, in particular in the response of the collective control inceptor. This inceptor, as one would  . 11 Phase diagrams of the collective and cyclic rotation response to heave acceleration input applied to the cockpit floor, at different excitation magnitude a, in m s −2 expect, is the most sensitive to vertical accelerations because it rotates in the vertical plane, acting on the rotorcraft heave axis control. When the complete spine model is considered, the resonance frequency of the collective response increases by 33%. A small rotation of the cyclic control inceptor, caused by the heave motion of the cockpit, has been observed. In both the longitudinal and lateral cases, the amplitude of the resonance peak is amplified when the dynamics of the spine is introduced; the resonance frequency increases by about 10% for the lateral cyclic response and by about 18% for the longitudinal cyclic.
The phase diagram of each lever's rotation angle is plotted in Fig. 11. The ratio between the instantaneous command rotation θ i and the deviation with respect to the trim configuration Δθ i is graphed against the analogous ratio on the command angular velocities. Both quantities are scaled with respect to the input acceleration level a. The results are obtained imposing a sinusoidal excitation to the cockpit, at 2.00 Hz, increasing the amplitude of the input vertical acceleration from 0.3 g to 0.6 g and 0.9 g. This frequency has been selected because it is very close to the resonance frequency for both the collective and the longitudinal cyclic dynamics.
A nonlinear dependence between the magnitude of the forcing term and the response can be observed in the plots of Fig. 11. This nonlinear behaviour is more pronounced in the cyclic response, whereas the collective only marginally departs from linearity. Figures 12 and 13 show several muscular activations on the left and right limb, respectively, imposing at the floor of the cockpit sinusoidal vertical accelerations at 2 Hz. From these graphs, it is evident that the maximum required activations increase with the amplitude of the forcing term. Moreover, as one would expect, the activations required for the left limb are much higher than for the right limb. For instance, the Abductur Pollicis Longus is the most active muscle in the left limb, with a maximum required activation of 0.26 at 0.9 g, whereas the maximum activation of the same muscle on the right limb is 0.03; an analogous consideration can be done for the Biceps Caput Longus. The Latissimus Dorsi shows a nonlinear behaviour as the amplitude increases: in particular for the left limb, a second harmonic component can be clearly observed for an excitation amplitude of 0.6 g. On the right limb, for an excitation amplitude of 0.9 g, the activations show sub-harmonic contributions at half the frequency of the excitation. This effect is most notable, for example, in the Deltoid anterior and in the Infraspinatus. For excitation amplitudes of 0.3 g and 0.6 g, the activations show a nearly sinusoidal response. This lower frequency component is related to a compound rotation of the cyclic lever about the longitudinal and lateral axes (Cf. Fig. 15) that is present at this specific frequency, while it is not noted at neighbouring ones: it is not exhibited, for example, at tests   conducted with an excitation frequency of 1.75 Hz. The phenomenon shows also a strong dependence on the backrest reclination angle: a moderate angle of reclination (about 10 • is typical in helicopters) will suppress it. It should be highlighted, in any case, that the test conditions were particularly harsh: a heave vibration, with magnitude close to 1 g, sustained for several seconds, can be definitely regarded as rare, in rotorcraft normal operation. Furthermore, no retaining devices, i.e. the so-called trim springs, were introduced on the flight control inceptors: the inceptors were left in the trim release configuration, in which only a balancing element, counteracting the effect of weight but offering no restoring moment, was active. Figure 14 shows some intervertebral moments about the medial axis between the two adjacent vertebrae for four periods of excitation at 2 Hz. As for the activation case, the moments increase as the forcing amplitude increase, but in particular for the moments between L5-S1 and T2-T3 the increment is nonlinear.

Conclusions
A multibody model of the upper body is proposed, obtained adding a detailed three-dimensional model of the spine to a previously formulated, detailed model of the upper limbs. A scaling procedure of the entire  Label K a (N m −1 )K bx (N m rad −1 )K by (N m rad −1 )K bz (N m rad −1 )C a (N s m −1 )C bx (N m s rad −1 )C by (N m s rad −1 )C bz (N m s rad − body model is to evaluate the biodynamics of helicopter pilots and occupants, for vibratory and rotorcraft-pilot couplings analysis. The proposed model shows good correlation with experimental data in terms seat-tohead transmissibility. Biodynamic feedthrough analysis of the complete upper body model shows that when the dynamics of the spine is considered, the frequency and the magnitude of the first resonance peak increase for all control inceptors; the largest increase, about 33%, is observed for the collective control incep-tor. Additionally, the cyclic control inceptor response departs from linearity when the amplitude of the acceleration imposed to the cockpit reaches 0.6-0.9 g at resonance. The capability to model the inverse kinematics and the inverse and direct dynamics of subjects performing various tasks, and the possibility to directly extract reduced order models to be coupled with vehicle models in linearised analyses makes the proposed analysis a powerful tool for comfort and pilot-vehicle interaction investigations.