The dynamics of human walking has been the subject of intensive studies since the ancient Greek period [1]. Mechanism models and parameters optimization are essential tools for the studies of human walking. These approaches are applied for a large variety of scenarios, including gait disorder simulation for a more and better understanding of performance and acceptable loads.
Multibody dynamics and methodologies are ideal for examining the biomechanics of human motion. The skeletal structure of humans comprises bones interconnected by physiological joints, which can be analogously represented as a conventional mechanism consisting of rigid links and mechanical joints. The human mechanism operates through muscle-tendon units, and it is governed by the central nervous system. Representing the human body as a mechanical multibody system enables the derivation of its equations of motion through conventional mechanical engineering methods. To accurately simulate human movements, these equations must be supplemented with models accounting for muscular activation and neural control [2].
Human walking has been analyzed employing different approaches based on optimization techniques to develop controllers for the simulation. A three-dimensional whole-body model was developed in [3] to predict human walking on level ground, and then used to explore the control strategies. Reinforcement learning was used to develop a torque-driven control for simulating walking and running [4], using also musculoskeletal models [5].
Muscle mechanical work in fact plays a crucial role in studying human movement: already existing experimental techniques struggle to precisely measure it. Utilizing muscle-driven forward dynamics simulations offers a reliable means to consistently analyze various aspects of mechanical work, from muscle fiber behavior to overall body movement. The approach proposed in [6] allows for an exploration of the demands of human locomotion, shedding light on both the complexity of muscle work and the constraints of traditional experimental methods.
In [7] a methodology was developed for generating muscle-actuated simulations of human walking for reproducing experimental measures of kinematics and ground reaction forces. In [8] the authors introduced a framework for simulating the interaction of muscle, ligament, and joint contact forces within the context of dynamic multi-joint movement. In a forward dynamic simulation, the differential equations of motion are numerically integrated forward in time subject to gravity, inertial and velocity-dependent effects, and muscle forces. Understanding the forces exerted on ligaments, tendons, and joints in typical conditions is vital for clinical purposes. However, acquiring such data often requires invasive procedures. Computer simulations offer a non-invasive manner to comprehensively monitor parameters like joint contact loads and soft tissue forces. Through simulation, the musculoskeletal system can be instrumented improving and moving beyond the feasibility with live human subjects [9].
Numerical techniques have been developed and proposed to solve problems related to many domains including engineering science, to which classical models are unable to provide reliable answers in reasonable times.
The term multi-Body generally refers to a set of bodies, rigid or deformable, connected via constraints, of which the relationship between acting/applied forces/torques and movements is studied in a wide range of conditions. Referring more specifically to the biomechanics of human body, the motion capture, is largely used to get the motion of the body segments that can be directly used to perform basic kinematics studies, [10–12]. The main goal is the creation of a realistic model for biomechanical analysis.
Presently, several software solutions exist, including both open-source and commercial options, offering such biomechanical models. Examples include in addition to the tools used in this paper, OpenSim [13] developed by NCSRR in Stanford, USA, BoB [14] created by BoB Biomechanics in Bromsgrove, UK, and AnyBody [15] from AnyBody Technology A/S in Aalborg, Denmark. These tools primarily function as post-processing tools, lacking the capability to deliver results in real-time.
The finite element method (FEM: Finite Element Method) is a numerical technique originally developed to solve problems related to industrial, civil and aeronautical problems in mechanics. The great potential of the mathematical model in solving complex non-linear systems in a large variety of conditions has meant that over the last 50 years the FEM method has become prevalent in the field of mechanical design, to study problems ranging from structural analyses, to dynamics, fluid dynamics, to name a few, while improving product quality and process and reducing production costs and times.
A field that has particularly benefited from these technologies is biomechanics, which refers to the discipline that uses the knowledge provided by mechanics in the study of the human body and its movements, and which finds its applications in various sectors:
-
In sport, through the study of the quality of the athletic gesture, it allows athletes to improve their performance while reducing the possibility of injuries.
-
In industry, to adapt work patterns and machinery for use by staff, assuming relevance in the innovative field of collaborative robotics (Cobot).
-
In medical science, such as orthopedic surgery and traumatology, in which the study of the locomotor system, both in healthy and traumatic conditions, associated with the study of the performance of prosthesis and implant solutions, is extremely useful option for the best possible treatment for the patient. Multibody biomechanical models can be also used in the forensic field to reconstruct and simulate tragic events like falls from a height, as described in [16]. In this work, a human body model made of 15 ellipsoidal elements, connected by 14 joints was developed in Adams.
In the last decades, the progress of medical science in the prosthetic field has been achieved also though the implementation of FEM method to study the response to loads and deformations applied to the complex geometry and structure of the bone, allowing the optimization of the design of prostheses and mechanisms fixation for fractures [ 17–19].
Various approaches are available in the literature for the modeling and simulation of complex systems, also involving the interaction of different fields, according to a multi-physics approach. Simulation tools give the possibility of analyzing gaits and postures as results of injuries in sports. In [20] the crouch walking with problem in lower limb shows high requirements of energy to make their body forward.
Pursuing competitive goals, closely related to innovation and technology, determines a compelling need to integrate digital processes into the design cycle of components, subsystems, and complex systems. This integration allows for comprehensive exploration of various aspects enhancing the correlation between physical and virtual tests, optimizing processes and facilitating efficient and effective product development in many sectors. Main and most famous example is co-simulation applied in autonomous driving as described in [21].
Co-simulation techniques begin to be implemented for a large variety of multi-physical problems. In [22] the design and simulation of an exoskeleton was carried out by developing a muscle-skeleton system designed in OpenSim, using mechanical parts of the exoskeleton developed in SolidWorks and then controlled by Matlab.
Co-simulation has been used in [23] to integrate exoskeletons and power tools into working environments, in which it is becoming frequent to use tools, like CoBots, automatic systems and operators sharing the same working area. In [23] OpenSim is used together with CAD software for the biomechanical analysis. In [24] a co-simulation framework was proposed including biomechanical human body models and wearable inertial sensor models to analyze gait events dynamically. Co-simulation techniques [25] were used for the analysis of grasping of deformable objects, but also for landing gear retraction and crash tube behavior due to hard and crash-landing analysis [26].
Experimental evaluation of human gait can be an efficient tool to fed into the models for getting realistic results of the motion or validating dynamic models, Data can be obtained by different technologies, such as electromyography (EMG), inertial measurement units (IMUs) and optical motion-capture systems.
Motion capture systems based on optical measurements are widely used for human gait analysis, efficiently used for kinematic analysis in real time, and they can be used for advanced analysis including torques by using exoskeletal models allowing the estimation of joint-reaction forces and muscle efforts. In [27] an extended Kalman Filter for optical motion capture system was proposed for estimating the kinematics and then with a multibody formulation evaluating the joint torques and muscle efforts producing such torques.
Wearable IMUs (inertial measurement units), alongside various measurement techniques, give the possibility to quantify human kinematics contributing to broadening the understanding of the complexity of the human movement in natural conditions [28] and in case of hemiparetic patients [29].
In this paper we use experimental data obtained by [30, 31] for the dynamic simulation of a human body model realized in MSC Adams environment and described and simulated in section 2, then the FEM of hip prosthesis is designed with MSC Marc and described in section 3. The introduction of co-simulation techniques is provided in section 4, while section 5 reports results of the numerical analysis.
Multibody modelling of the human walking
The extent of detail in a multibody dynamic model for human movement hinges on the specific predictions and analytical outcomes sought. Human motion can be simulated using ideal joint torques, activations, or muscle excitations. Some predictive simulation studies have used conceptual models with minimal degrees of freedom (DoF) and without muscle representations. Alternatively, musculoskeletal (MSK) models or more complex neuromusculoskeletal (NMSK) models, incorporating the central nervous system (CNS), may be employed to capture additional dynamics. In the following skeletal representation is considered.
Multibody models depict the movement of individual segments of the human body through rigid or flexible bodies. Each segment encapsulates a bone or a cluster of bones along with the encompassing tissue that moves accordingly and coherently. The connections and interactions among segments are characterized by constitutive laws that emulate elastic components like joint cartilage and ligaments. Kinematic constraints may also be applied to articulate the relative motion between segments.
The skeletal model can be activated using an ideal joint torque. When employing an NMSK model, it can be activated through muscle force, activation level, or muscle excitation controlled by a neural controller.
The typical representation of the kinematic and dynamic model of the human body involves a sequence of rigid segments or articulated linkages interconnected via joints [32]. From a system topology perspective, the skeletal system primarily consists of an open kinematic chain, although less frequently, it could define a closed kinematic loop (e.g., when both hands grasp a rigid object). An open kinematic chain model has been employed to simulate the swing phase of human walking [33]. The open-loop model, characterized by serial kinematic connections and employing a minimal number of coordinates, gives rise to a system of ordinary differential equations (ODEs) [34]. Conversely, when utilizing dependent coordinates (e.g., natural coordinates), it yields a system of differential algebraic equations (DAEs).
The Body segment configurations are delineated by generalized coordinates, which may take on either relative or absolute attributes. Relative coordinates chiefly pertain to internal motion, depicting the relative displacement between successive segments, whereas absolute coordinates signify position and orientation concerning an inertial reference frame. These generalized coordinates serve as the foundation for defining the system states. Bodies are typically presumed to be rigid, with segment properties such as length, center of mass (CoM), mass, and moment of inertia often derived from anthropometric measurements. These measurements, sourced from anatomical regression equations found in the literature, are frequently adjusted to correspond to the specific characteristics of the subject under study. experimental data from motion capture systems may be used to estimate segment lengths. For fast motions, the segments could be modeled as wobbling masses [35].
In human movement modeling and analysis, ideal joints serve to link adjacent rigid segments, introducing degrees of freedom (DoF) and thereby imposing constraints on the system. These joints can range from simple 1-DoF revolute (pin) joints [36] to more complex 6-DoF joints [37]. Kinematic constraints play a crucial role in modeling anatomical features that are challenging to represent using standard joints or would necessitate intricate ligament and cartilage models to accurately capture motion characteristics [38]. Models of anatomically complex joints such as the knee, shoulder, and neck often rely on kinematic constraints to achieve the desired level of detail [39].
The prevalent external contact typically addressed and incorporated into multibody human models is the interaction between the foot and the ground. This aspect is essential for accurately simulating activities such as walking, running, and jumping. Contact can be modeled in a variety of ways. A rigid contact model might overlook certain foot movements occurring during contact. One enhancement involves incorporating the "roll" motion of the foot, achieved by modeling the foot as a rigid body capable of rolling without slipping. Apart from rigid contact, various deformable contact models exist. Among these, point-contact models stand out as the most straightforward. In a point contact scenario, the normal forces are typically represented by linear or non-linear springs and dampers connecting the ground to a single point on the contact surface. An alternative approach, known as volumetric contact, employs analytical geometry to characterize surfaces. This model utilizes elastic foundation theory to distribute contact loads across a surface, rather than confining them to discrete points. Initial explorations of volumetric contact in foot-ground interactions have employed shapes such as spheres and ellipsoids, while its application has extended to modeling impacts in various sports contexts [40].
Segment properties such as mass, length, and inertia, along with joint centers and axes, are typically customized from a generic model to a specific subject model using regression equations. This approach is common in most predictive simulation studies [41], where segment lengths are measured using reflective markers and utilized to calculate scaling factors for adapting anthropometric properties. Depending on the study's objectives, population-specific models could be employed, such as those tailored for different age groups or genders, like young adults versus older adults, or females versus males. However, it's worth noting that linear scaling methods have certain limitations and may not be accurate enough when subjects are far away from the mean population values.
Some studies utilize dynamic measurements, such as marker trajectories and ground reaction forces, in conjunction with optimization techniques to refine model joint and/or inertial parameters. For instance, in [42], a methodology was devised to estimate hip, knee, and ankle joint centers and axes, along with full-body inertial parameters. Building upon this approach, the authors predicted the walking patterns of individuals with osteoarthritic knees using a 3D model driven by ideal torques, aiming to tailor gait modifications specific to knee osteoarthritis rehabilitation [43].
Taking into consideration the limitations of linear scaling methods, in this paper we have considered a multibody model using MSC Adams by Hexagon realized according to the experimental data available in [30–32]. This model is composed by rigid segments, each with proper mass and inertial properties, and connected through elastic joints. The initial conditions for each part of the human body model, in terms of position and velocity of the centers of mass, were determined through a design of experiments (DOE). The analysis was based on the final body position, as depicted in pictures and measurements provided by legal prosecutors. [31–32].
Experimental data regarding the walking motion of the human body was provided by [30]. The procedure describing the implementation of the biomechanical model is available in the work from [27]. The test environment consisted of a typical Motion Capture setup, with 36 optical markers fixed on a lean and muscular-built human male, and two force plates on the floor measuring the foot-ground contact forces, as reported in [30]. The data includes the kinematics and topology data for a two-step, 1.25s long walk, and was implemented in a Matlab file for ease-of-use [30–32].
The biomechanical model consists of 18 rigid bodies, connected through 17 ideal spherical joints, thus defining 57 degrees of freedom in a Cartesian space, as depicted in Fig. 1. The computational model was defined with 228 mixed (natural + angular) coordinates, including, for each time increment the natural coordinates of joint centers, optical markers and the body segments’ centers of mass. In Fig. 2, the hollow white dots represent the filtered optical marker positions, the full gray dots are the obtained joint centers while the triaxial markers represent the COM relative coordinate system (red arrow for x axis, green arrow for y axis, blue arrow for z axis). The latter were obtained from the rotation matrices and the joint centers positions, according to the following
$$\:\overrightarrow{{P}_{G}}=\overrightarrow{{P}_{J}}+\stackrel{̿}{R}\:\overrightarrow{d}$$
1
Where:
-
\(\:\overrightarrow{{P}_{G}}\) represents the vector containing the coordinates for the center of mass, relative to the current body segment and time increment, in the absolute coordinate space.
-
\(\:\overrightarrow{{P}_{J}}\) stands for the vector containing the coordinates for the proximal joint center, relative to the current body segment and time increment, in the absolute coordinate space.
-
\(\:\stackrel{̿}{R}\) is the rotation matrix
-
\(\:\overrightarrow{d}\) is the distance vector between the joint center and the center of mass in the initial, t = 0s, configuration.
The Matlab model was used for a graphical representation of the motion, including the force vectors on the plates [30–32], having 0.01s time increments.
The multibody software used in this work is MSC Adams 2022.3. The Adams software uses the Lagrangian (aka energetic) approach to the solution of dynamics PDE.
The biomechanical rigid-body model in MSC Adams was built by referring to the initial configuration reported in Fig. 2. To consider the difference in the reference systems, a transformation matrix was used. Particularly, in the Adams model X axis refers to the antero-posterior direction; the Y axis refers to the vertical direction, and the Z axis refers to the medial-lateral direction, with opposite sign relative to the Matlab model’s Y axis. Overall, 17 points were defined in the t = 0s positions of the joint centers, then 18 rigid bodies of suitable characteristics were modeled. The bodies were then connected by 17 ideal spherical joints. This is good for a first approximation, as real joints do not possess a full rotation about all axes, and do not act as perfect locks, but react to the relative movement with a strongly non-linear, viscoelastic behavior. For the knee, although the center of the joint noticeably moves during the flexion-extension, for the purpose of the current work, it was decided to neglect this phenomenon. The foot was modeled as made of two pieces, connected by a spherical joint at the metatarsal, to improve the foot-ground contact force. The mass and inertia information about the body segments were obtained by [Coruna], which portraits the human body mass properties of a 1,75m tall human male weighing 75kg and reported in Fig. 3. The inertia tensor components are considered relative to the body segment’s center of mass.
A topological map of the model is reported in Table 1, while body parts characteristics are given in Table 2.
Table 1
Multibody human body joints
| JOINT | DESCRIPTION |
|---|
| JOINT_1 | Lumbar joint |
| JOINT_2 | Right hip joint |
| JOINT_3 | Left hip joint |
| JOINT_4 | Neck joint |
| JOINT_5 | Right shoulder joint |
| JOINT_6 | Left shoulder joint |
| JOINT_7 | Cervical joint |
| JOINT_8 | Right elbow joint |
| JOINT_9 | Right hand joint |
| JOINT_10 | Left elbow joint |
| JOINT_11 | Left hand joint |
| JOINT_12 | Right knee joint |
| JOINT_13 | Right ankle joint |
| JOINT_14 | Left knee joint |
| JOINT_15 | Left ankle joint |
| JOINT_16 | Right metatarsal joint |
| JOINT_17 | Left metatarsal joint |
Table 2
Multibody parts characteristics
| BODY | Type | MASS [kg] | Ixx [kg/mm^2] | Iyy [kg/mm^2] | Izz [kg/mm^2] |
|---|
| Pelvis | Plate | 11.28 | 3.3E + 05 | 4.9E + 05 | 3.3E + 05 |
| Trunk | Plate | 27.36 | 1.22E + 06 | 8.6E + 05 | 8.9E + 05 |
| Neck | Link | 1.7325 | 1.01E + 04 | 1.01E + 04 | 1.01E + 04 |
| Head | Sphere | 3.465 | 6.0E + 04 | 3.22E + 04 | 6.0E + 04 |
| Right_arm | Link | 2.415 | 6.8E + 04 | 8500.0 | 6.8E + 04 |
| Right_forearm | Link | 1.32 | 1.75E + 04 | 4700.0 | 1.75E + 04 |
| Right_hand | Link | 0.4875 | 4000 | 1200 | 4200 |
| Left_arm | Link | 2.415 | 6.8E + 04 | 8500.0 | 6.8E + 04 |
| Left_forearm | Link | 1.32 | 1.75E + 04 | 4700.0 | 1.75E + 04 |
| Left_hand | Link | 0.4875 | 4000 | 1200 | 4200 |
| Right_leg | Link | 7.11 | 2.8E + 05 | 6.7E + 04 | 2.8E + 05 |
| Right_shank | Link | 3.2025 | 1.34E + 05 | 1.81E + 04 | 1.34E + 05 |
| Right_foot | Plate | 0.9 | 193.53380899 | 192.32625613 | 2.0332932482 |
| Left_leg | Link | 7.11 | 2.8E + 05 | 6.7E + 04 | 2.8E + 05 |
| Left_shank | Link | 3.2025 | 1.34E + 05 | 1.81E + 04 | 1.34E + 05 |
| Left_foot | Plate | 0.9 | 202.00239037 | 200.74199785 | 2.1222653504 |
| Right_toe | Link | 0.105 | 3.190779708 | 3.1708708601 | 3.35227775E-02 |
| Left_toe | Link | 0.105 | 3.1962505943 | 3.1763076109 | 3.35802555E-02 |
For simulation purposes, a floor, consisting of a 5x0.02x5m element was designed and the contact forces with the foot segments were defined using a Hertz normal impact and Coulomb friction formulation.
MSC Adams allows to create joint motions for the ideal spherical joints. These motions are defined with a Body123 method, namely the first rotation is performed according to the body’s X axis, the following is performed according to the Y axis of the new, rotated, relative coordinate system, and the last rotation is performed according to the Z axis of the previously twice rotated relative coordinate system. The rotation matrices, provided by the Matlab file, are relative to the centers of mass in a standard, space fixed, formulation, and must be transformed to define joint motions in the Adams model. Instead, it was decided to add 17 general motions, between each body segment and the ground, located in the joint centers, prescribing all coordinates. The joint center position history was imported, from the relative Matlab matrix, into the Adams model, in the form of time-based splines, one for each component, through the Import test data interface window. The motions were then defined by prescribing the translational coordinates of each joint, avoiding to over-constrain the model, and using the AKISPL formulation. The motion is reported in Fig. 3. Since prescribing translational motion for a rather complex multibody system typically makes it hard for the Solver to converge, the ideal joints were suppressed for this analysis, and substituted by bushings, which then acts as real joints, with compliance, allowing small relative motion.
A dynamic analysis, consisting of 620 integration steps for a total simulation time of 1.24s, was then performed. Applying translational motion is, in fact, an approximation, as the vertical translations are constrained, and the model is not moving because of muscle activation and the friction between the feet and floor, but due to a given motion. This analysis was then only carried out to obtain the rotations of each body segment, with the appropriate body123 formulation, through orientation measures, in a dynamic environment. The ideal joints were then re-activated, all but two, and provided with the relative joint motions, using the splines obtained from the orientation measures (Fig. 4).
To avoid an over-constrained model, while giving it some degree of self-adjusting dynamic capabilities, the ideal joints, and motions, relative to the metatarsals (JOINT_16 and 17) were suppressed and substituted with bushings, then providing a passive behavior closer to that of a real joint and allowing the toes to better adhere to the floor without expressing excessive contact force peaks. The orientation measures could not provide data relative to the most peripheric bodies, including the head and both hands. The relative motions were then suppressed, imposing a 0*time displacement: the observed small angular displacements and inertia tensor components of the parts made them considered as negligible. The three plots in Fig. 4 represent, respectively, the x, y and z components of the translational displacement of each body part’s COM, obtained in the translational motion-driven simulation. In Fig. 5 the nomenclature used for the rotational measures is explained as follows referring to the JOINT_xx, with xx being the joint number as noted in Table 1 and Fig. 3; the Bxx_1, Bxx_2 and Bxx_3 represent, respectively, the first rotation about x; the second rotation about the new, rotated axis y’; the third rotation about the new, rotated twice axis z’’.
To provide some qualitative considerations, as expected, the translational displacement components are prevalent on the x-axis of the sagittal plane, along which the movement is performed, while on the y and z axes the range of movement is smaller, as the vertical displacement is a result mainly of the flexion-extension of the leg segments, and the transverse one is mainly due to the arms and legs’ adduction-abduction. As for rotational displacements, the higher values are registered around the z-axis, representing the flexion-extensions, while the range for the motion around the x-axis, standing for the adduction-abductions, and the y-axis, representing rotations, are predictably lower.
While the filtering formulation used to obtain the joint positions in [30] may have suppressed this to a certain degree, it is to be noted how the right and left side of the model exhibit a slight asymmetry: since the data are experimental acquisition-based, this is to be expected, as the human body is naturally asymmetric, even for healthy subjects.
A set of 4 three-component bushing-like torques was added, to the trunk and pelvis body segments, with the objective to act as a control system for the body. The human body has a complex neuro-muscular control system that uses multiple sensory inputs to maintain balance and move properly. In the Adams dynamic environment that was defined here, some sort of simple control system was also required for a successful simulation.
Two gain parameters and angular displacement and velocity error measures were created. These were defined as the difference between the current angular displacement/velocity measures and the splines obtained from the corresponding measures in the posgiunti4 analysis. The torque components are then defined in order to act as a bushing-like PID system. It was necessary to gradually increase the nominal value for the torque components in the initial steps of the analysis to avoid sudden spikes, so, a quintic polynomial approximation of the Heaviside step function was used to reach the nominal value of the torque in a small but finite amount of time. The function syntax is reported in (2).
| \(\:\varvec{i}\in\:[\varvec{x},\varvec{y},\varvec{z}]\) | Referenced axis of the coordinate system |
|---|
| TQ1i(t) | Components of the torque based on angular displa cement error |
| TQ2i(t) | Components of the torque based on angular velocity error |
| G1,G2 | Gain parameters for the two torque sets |
| t1,t2 | Time values at which the step function ends |
| t10, t20 | Time values at which the step function begins |
| EAi(t) | Error components based on the angular displacement |
| Eωi(t) | Error components based on the angular velocity |
$$\:\varvec{\alpha\:}{\varvec{n}}_{\varvec{i}}\left(\varvec{t}\right)=\left\{\begin{array}{c}Gn*E{A}_{i}\left(t\right)-h{n}_{0}\:\:\:\:\:\::n=1\\\:Gn*E{\omega\:}_{i}\left(t\right)-h{n}_{0}\:\:\:\:\:\::n=2\end{array}\right.\:\:\:\:\:\:;\:\:\:\:\varvec{h}{\varvec{n}}_{\varvec{i}}\left(\varvec{t}\right)=\left\{\begin{array}{c}Gn*E{A}_{i}\left(t\right)\:\:\:\:\:\::n=1\\\:Gn*E{\omega\:}_{i}\left(t\right)\:\:\:\:\:\::n=2\end{array}\right.$$
$$\:\varvec{\Delta\:}\varvec{n}\left(\varvec{t}\right)=\frac{t-t{n}_{0}}{tn-t{n}_{0}}\:;n\in\:\left[\text{1,2}\right]$$
$$\:\varvec{T}\varvec{Q}{\varvec{n}}_{\varvec{i}}\left(\varvec{t}\right)=\left\{\begin{array}{c}h{n}_{0}\:\:\:\:\:\:\:\::t\le\:t{n}_{0}\\\:h{n}_{0}+\alpha\:*{{\Delta\:}\text{n}}^{3}\left[10-15{\Delta\:}\text{n}\left(\text{t}\right)+6{{\Delta\:}\text{n}\left(\text{t}\right)}^{2}\right]\:\:\:\:\:\:\:\::t{n}_{0}<t<tn\\\:h{n}_{i}\left(t\right)\:\:\:\:\:\::t>tn\end{array}\right.\:\:\:;\:\:\:n\in\:\left[\text{1,2}\right]$$
$$\:T{Q}_{i}\left(t\right)=\sum\:_{n=1}^{2}TQ{n}_{i}\left(t\right)$$
2
A new set of dynamic analysis (vforce4_2 in results) was then performed. The contact impact and friction, and PID torques gain parameters, were subject to a series of numerical optimizations, with the objective of improving the dynamic behavior of the model and its coherence with the experimental source kinematics. The consistency of the behavior of the model was then verified.
The following results were obtained, as reported in Table 3 and Fig. 7.
Table 3
| Stiffness (N/mm) | Exponent | Damping (Ns/mm) | Max penetration (mm) |
|---|
| 300 | 1.0 | 6 | 0.1 |
| µs | µd | Stiction transition velocity (mm/s) | Friction transition velocity (mm/s) |
|---|
| 0.7 | 0.5 | 232.9 | 2329.0 |
A comparison of the Body123 body rotation measures, between the analysis with the translational motions and the rotational motions one shows that, since the motion on both feet and hands was suppressed, and their rotation is only determined by the interaction between the contact forces and the bushing joints, the only noticeable differences are on the head joint, and on the left knee and ankle.
By comparing the results, the model kinematic behavior can be considered consistent. Particularly, in Fig. 8, the vertical component from the total contact force on all foot segments (CONTACT_T_y) is compared to the vertical components measured on the force plates obtained as experimental result of the MoCap (FUNCTION_MEA_9 and 10). It is worth noting that the regime values from the simulation are consistent with the measured data, except for peak values and transients, which are the consequence of the default impact contact formulation. The Hertz impact function is valid for materials exhibiting a linear-elastic behavior, while the human foot (eventually wearing shoes) has strongly non-linear characteristics and is then not properly represented by the rigid plate body used in this context. Further studies on this subject are to be carried, but for the purposes of this study the assumption can be considered acceptable, since the main goal is not modelling the contact with the ground or impact of the foot.
The rotational displacement total error, between the translational displacement motions-driven analysis and the new ones, was also provided.
The plot in Fig. 9 shows a comparison between the rotational error of an analysis where only 1 torque PID was provided on the trunk (red full line), and the final one where there are 4 torque PID’s (blue dash line).
The benefit is evident, with the maximum total error dropping from 0.8° to 0.65°, and the average and final behaviors improving greatly. It appears even more so if considering the total torque applied to the model by the PID system, as reported in Fig. 10.
The peak values, at about 10^6 Nmm, are comparable, while the 4 torque PID analysis has a better behavior.
The Adams multi-body model was then deemed as acceptable, accepting the validation using the experimental data provided and for the purposes of this study, and was then prepared for the MBD-FEM co-simulation.
FEM model synthesis
FEM analysis techniques are being used in the biomechanics field, since the late 70’s, to provide a feasible method for the study of the human tissue, to improve both processes and products.
In the case of joint replacement operations, the general trend for an increased life expectancy determined the need for better performing prosthesis, which can last more, with less aftereffects, and require less invasive maintenance operations. In this sense, it is desirable to design a method that, implementing several technologies, goes from the medical screening of the particular patient to the ad-hoc designed and optimized final product.
In this study, a hip prosthesis designed and reported in [44], was inserted into a 3D human femur model (Figs. 11 and 12), imported into MSC Marc, meshed and subject to the displacement field provided by the CoSim interface. In this way, instead of using the static tests prescribed by the current standards, the loads experienced by the FEM model can be designed as more similar real dynamic conditions, coupled with a gait analysis.
In addition, such procedure can provide a prosthesis that is tailored to the specific patient’s needs.
The Marc .dat model consists of two meshed solids, one representing the prosthesis, the other representing an average human male femur, subjected to Boolean cut operations to provide the space for the prosthesis.
Given the quite complex and strongly asymmetrical nature of the geometry and boundary conditions, a 3D structural analysis needs to be performed. The solids were meshed through the Automesh volumes Marc function, resulting in two sets of tetrahedral full integration elements.
The geometric properties for all elements included the Assumed strain option to improve shear behavior. No part was supposed to act in plastic regime, so it was not necessary to flag Constant dilatation. Mesh data is given in Table 4.
Table 4
| Element set | # of elements |
|---|
| protesi (prosthesis) | 9442 |
| femore (femur) | 31796 |
| Nodes set | # of nodes |
|---|
| protesi_n | 2498 |
| femore_n | 7110 |
The corresponding material properties were defined as: for the prosthesis, the linear elastic data for the Ti6Al4V alloy were provided. It was decided to model a linear elastic behavior to the femur too, ignoring the actual viscoelastic behavior, and the differences between cortical and trabecular bone, according to the limitations imposed by the quasi-static analysis necessary for the CoSim interface (Table 5).
Table 5
| Material | Young modulus [GPa] | Poisson ratio |
|---|
| Ti6Al4V | 113 | 0.31 |
| Bone | 20 | 0.3 |
Two contact bodies were also defined, interacting through an ideal glue formula, the strongly non-linear behavior of the cement, used in actual hip replacements, could not be provided due to lack of experimental data. Instead, it was decided to impose a glue, rigidly fixing the two contact bodies together.
It is to be noted that Marc does not require the user to specify the units for the various physical quantities: only the length unit is requested, to properly scale CAD imports. The a-dimensional physical quantities are then used in the calculations: for example, stress is expressed as “M L^-1 T^-2”. It is standard procedure, then, to choose the most appropriate units while importing data, according to the desired results: with the main goal to obtain MPa for the stresses, while having mm as length unit, the t (metric ton) mass unit is usually considered.
For a co-simulation analysis, the CoSimPRE interface requires the user to specify the Marc model units to ensure proper communication with the Adams model. It was found that it is necessary to consider the same MMKS (mm, kilogram, second) unit system, used in Adams, also in Marc. The stress unit expressed in the result files is then to be considered as \(\:\frac{kg}{mm*{s}^{2}}=kPa\).