FE Model and Operational Modal Analysis of Lower Limbs

Human lower limbs are exposed to vibrations on a daily basis, during work, transport or sports. However, most of the FE (Finite Elements) and OMA (Operational Modal Analysis) studies focus either on the whole body or on the hand-arm system. The study presented herein aims at identifying the modal parameters of the lower limbs using a 2D FE model updated using OMA. A numerical model is proposed, and a modal analysis has been performed on 11 subjects. Two repeatable modal frequencies were extracted: 52.54± 2.05 Hz and 118.94± 2.70 Hz, which were used to update the mechanical properties of the numerical model. The knowledge of these modal characteristics makes it possible to design new equipment that would absorb these specific vibrations and possibly reduce the risk of related diseases in the field of sports and transport.


Introduction
Every day, the human body is exposed to vibrations, whether during work, transport or sports.Two modes of exposition can be distinguished: whole body vibration and local vibration.Whole body vibrations may have some beneficial effects including ossification [1,2] and treatment of child muscular disabilities [3].However, their negative effects have also been proven.Vibration could lead to an accelerated degeneration of joints [4] and vascular diseases [5].Likewise, local vibrations produce beneficial effects, in sports training for example [6], as well as significant drawbacks, especially in the hand-arm system.Related symptoms can be split into four groups: vascular effects, muscular effects, skeletal effects and joints effects [7], including Raynaud's disease [8].All of them are generally referred to as Hand Arm Vibration Syndrome (HAVS) [9].
Investigating the effects of vibrations leads us to consider the human body as an extremely complex system, if only because of its responses to vibrations.It is a non-uniform, nonlinear and anisotropic structure.Nevertheless, simplified numerical models, like Finite Elements (FE) models, have been used to study the human body on various aspects: evaluate the injury risk of the lumbar vertebra of a seated human body submitted to vibrations [10] or predict the strains and strength of human femora [11].A recent study featuring a two-dimensional FE model of the human hand-arm system was also used to estimate significant modal frequencies [12].These frequencies were then compared to experimental results in order to update the mechanical properties of the model.The study concluded that FE models could be used to obtain reliable resonance frequencies of the human hand-arm system.
In order to update and validate numerical models, it is necessary to carry out experimental measurements.In this study, Operational Modal Analysis (OMA) is used.The OMA represents the set of modal identification methods based on the sole measurement of the system's response.It has been successfully used in various studies on the human body, mostly for whole body vibrations.Indeed, one of the main fields of interest of modal analysis in this scope is related to transport, including the effects of cabin vibrations [13] or the resonance behavior of the seated human body [14].OMA and finite elements have already been used jointly to determine modal frequencies, damping and mode shapes, but only on a given part of the human body, e.g., the tibia [15].Moreover, most of the studies about vibrations and the human body only focus on the hand-arm system: the upper limbs are indeed frequently exposed, for example during biking [16].However, lower limbs should be given proper attention, as well, since many everyday activities expose legs to vibrations.Solicitations could be similar to white noise in public transport [17] or similar to shock in walking and running activities [18,19].Therefore, we have chosen to focus our study on the lower limbs, combining OMA and finite elements on a different scale, i.e., the half-body.
The aim of this study is to identify the experimental modal parameters of the human lower limbs and to propose an updated two-dimensional finite element model.Experimental modal parameters were extracted by operational modal analysis using poly-reference implementation of a Least Squares Complex Frequency domain (p-LSCF) method.

Measurement
Eleven voluntary and healthy subjects were recruited for this study.Their anthropometric parameters are summarized in Table 1: age, height, weight, Body Mass Index (BMI) and characteristic lengths.Height was measured with a stadiometer and weight with an electronic load cell scale, and BMI was calculated as the weight divided by the square of the height.All characteristic lengths were measured with a thin flexible metal tape measure.The study was approved by the local university ethics review board and in agreement with the Declaration of Helsinki.The participants were aware of the purpose of the study, and they all provided written informed consent.The exclusion criteria included a history of back pain, acute inflammation in the pelvis and/or lower extremities, acute thrombosis, recent fractures, recent implants, gallstones, kidney or bladder stones, any spine disease, peripheral vascular disease and severe delayed onset of muscle soreness in leg muscles.
The study was divided into two protocols.The first protocol (Protocol 1) was dedicated to the repeatability of the OMA results.It consisted of 10 repeated measures on the same subject (healthy man, age 20, 78 kg, 180 cm).The second protocol (Protocol 2) identified the modal characteristics (modal frequencies, damping and mode shapes) of the lower limbs of all subjects.Ten measurements were performed on each one of them.During the protocol, the subject had to be barefoot to create a single impact, making it possible to compute all of the resonant frequencies.To create the impact, each subject had to hold a rod and drop it from a height of 5 cm.Subjects were asked to keep their legs straight, to not damp the shock by bending their ankles or knees and to hit the ground using their heels.Subjects were finally instructed to limit muscle contraction.For the first protocol, the subject performed ten drops, and each of the measurements was independently analyzed for repeatability purposes.For the second protocol, eleven subjects performed ten drops each.These ten drops were considered as a unique set of data and analyzed accordingly.
Both protocols relied on the same experimental device, which was based on two synced Oros OR 35 multi-channel analyzers (Oros, Grenoble, France).The recording lasted one second, with a sampling rate of 2560 Hz.The impact response of the lower limbs was measured using four tri-axial accelerometers (Bruel and Kjaer, Naerum, Denmark) (Deltatron 4524 B; frequency range: 0.25-3000 Hz; sensitivity: 10 mV/g).Those were tightly fastened close to the joints in the lateral side of the right limb to minimize the contributions of skin artifacts, at repeatable anatomic locations: at the top point of the iliac crest, at the center of rotation of the knee placed in the middle of the lateral condyle, at the middle of the ankle external malleolus and at the center of mass of the tibia, aligned with the other two accelerometers on the segment.The accelerometers' axes were properly oriented in the same direction (Figure 1).Accelerations along the x, y and z axis were recorded at all locations.The distance between the accelerometers was measured in order to rebuild a geometry using Oros Modal Analysis software (Version 5, OROS, Grenoble, France).The distances between the hip and knee accelerometers, between the knee and tibia sensors and between the tibia and ankle sensors are summarized in Table 1.

Modal Parameters
Modal identification was performed using OROS Modal software (Version 5, OROS, Grenoble, France).We chose the Least Squares Complex Frequency domain (p-LSCF) method [20] owing to its capability of separating strongly-coupled modes and avoiding the decomposition of model-related residues into singular values.Singular Value Decomposition (SVD) drastically decreases the quality of the regression [21].The method is based on the analysis of the outputs, described in Equation (1): ω being the frequencies, H(ω) the matrix composed of the frequency response function, exponent H the complex conjugate transpose of a matrix (Hermitian) and S inputs and S outputs the inputs' and outputs' spectral matrix.The outputs' estimation S+ outputs (ω) is obtained by integrating two residual terms to approximate the effects outside the selected range [22], using Equation (2): where T is the transpose of a matrix, * is the complex conjugate of a matrix, Φ k is a column vector containing mode shape k, g T k is a line vector with the modal participation factor of mode k and λ k are the structure poles.LR and UR stand for the lower and upper residues, which illustrate the influence of modes outside the frequency range.
For both protocols, the range of interest was 0-200 Hz.Preliminary measurements showed that for frequencies over 200 Hz, the spectral signal was similar to noise.Moreover Standard ISO 2631-1 [23] provides a guideline for measurement and assessment in the case of exposure to whole body vibration, and defines a frequency weighting filter.More precisely, frequencies above 200 Hz are significantly filtered.
Identification of the modal parameters, and particularly mode shapes, requires a simplified model.Thus, a model of a segmented straight line was made from four nodes with three degrees of freedom (DOFs), positioned in a way similar to the measurement points in OROS Modal software (see Figure 1).For Protocol 1, the modal characteristics of the leg were calculated: modal frequencies, damping and shapes.The results of frequencies and damping are presented as the mean and standard deviation (coefficient of variation).The coefficient of variation was calculated as the ratio between the standard deviation and the mean.The repeatability of the shapes was visually estimated.For Protocol 2, the modal characteristics of the lower limb were determined for all of the subjects: modal frequencies, damping and shapes, as indicated above.The modal frequencies were then used to validate and update the mechanical properties (Young's moduli and mass densities) used in the numerical model.

FE Model
The 2D FE model of the human lower limb was created using Abaqus R software (Version 6.12, Dassault Systèmes, Vélizy-Villacoublay, France) and features the foot (calcaneus, talus, metatarsus lumped together, tarsus lumped together), leg (tibia bone), thigh (femur bone), soft tissues, ankle, knee, hip joints and trunk (Figure 2).The femur and tibia bones were made of a central part of cortical bone and trabecular bones at the ends.The bones were held together by soft tissues, which represent muscles, skin, tendons and adipose tissue.Bones are in contact at joints; however, joint cartilage has not been modeled.The anthropometric dimensions used to create the model have been retrieved from several references (Table 2) and calculated for a virtual subject 175 cm tall.Obviously, the natural modes and frequencies are strongly dependent on these anthropometric dimensions, yet focusing on a single virtual subject showed promising results and allowed us to validate the methodology.
The model was developed for a standing posture without knee flexion.Humans are indeed exposed to vibrations and impacts while extending legs.In industry, workers are liable to stand on platforms, which are sources of vibrations, and it is common to remain standing in public transport.
Four materials were used: trabecular and cortical bones and soft tissues for the trunk.All of the components of the trunk (spine, scapular, abdomen, etc.) were lumped together (see Figure 2) to simplify the model, as defined by Adewusi [12].Mechanical properties of the bones [24], soft tissues and trunk [12] are summarized in Table 3.The downside regarding the use of a linear elastic constitutive law for soft tissues is discussed in Section 4. Fixed boundary conditions were applied at the shoulder to allow motion of the trunk and hence the pelvis.The mesh is created using default CPS33-node linear plane stress triangle elements, for a total of 1167 nodes and 1742 elements.Modal analyses of the model were then performed using Abaqus R on a standard workstation.Modal frequencies' sensitivity to the mesh density was investigated.

Segment Parameters
Dimension (mm)

Modal Identification
The repeatability results are listed in Table 4. Eight modes have been identified in the selected frequency range.The score stands for the number of appearances of the mode compared to the total number of measures.Modes 1, 3, 4, 6, 7 and 8 have been identified in more than 80% of the tests.Mode 2 was identified in 60% of the tests and Mode 5 in 30%.The mode frequencies exhibit a coefficient of variation below 5%.Frequencies 1, 3, 4, 6, 7 and 8 have been identified in at least 80% of the tests, and their coefficients of variation lie below 5%.Their identification is therefore considered repeatable.Dampings for Modes 2, 3, 4, 6 and 7 have similar values, around 7.57%.Damping of Mode 1 is 40% higher, while that of Mode 5 is about 50% lower and that of Mode 8 only 30%.The variation coefficients range from 16.2% for Mode 6 to 49.9% for Mode 5.The modal dampings are consequently considered not repeatable, contrary to the frequencies.The parameters considered by the protocol could not provide accurate measurement of damping.Therefore, damping will be discarded in the forthcoming exploitation of the results.
The results of the modal identification for the 11 subjects are listed in Table 5.Eight modes have been identified among all of the subjects.The modal frequencies were estimated as identical for the eight identified modes with a coefficient of variation below 5%.However, modal damping had a high coefficient of variation (above 30.1%).This confirmed the previous results of the repeatability protocol.

FE Model Results
The results of the initial numerical model highlight three modes: 12.8 Hz, 54.5 Hz and 126.9 Hz, depicted in Figure 3.The deformations of the first mode are located at the top of the trunk and the foot.Experimental mode shapes have been visually compared to the numerical ones: the modes corresponding to the numerical modes were Numbers 3 and 8. Mode 3 matched numerical Mode 2; Mode 8 matched numerical Mode 3. The frequency of Mode 3 is 52.54 Hz and that of Mode 8 is 118.94Hz.Given that our study was based on a two-dimensional model, the other modes could not be obtained by the numerical model.The corresponding experimental mode shapes showed however that deformation is mostly oriented out-of-plane, suggesting that the numerical model accurately describes the in-plane dynamic behavior.

Model Update
The mechanical properties used during the modal analysis are summarized in Table 3.A sensitivity analysis was performed in order to investigate the effect of changes in the mechanical properties on the resonant frequencies [34].It was observed that the mechanical properties of the soft tissues were the most influential on the resonant frequencies of the FE models (a variation of 50% of the parameters of soft tissues, including the trunk and hard tissues, leads to a variation of about 10% in the frequencies for soft tissues and 1% for hard tissues).Then, the mechanical properties of the bones were retained while those of the soft tissues and the trunk were updated.Actually, the properties of the trunk were not known, and those of soft tissues depended on the subject stature and muscular force.The aim was to update the mechanical properties of the model until its resonant frequencies were close to the OMA results.The different frequencies obtained from the numerical model and with the OMA were associated through visual comparison of the mode shapes.
A Shapiro-Wilk test showed that the values of experimental frequencies for Modes 3 and 8 were normally distributed.A Student test has been used to calculate a confidence interval for each frequency, as shown in Equation ( 3 where t α was set up for a risk probability α of 5%, n was the number of subjects and fi exp and s i were the mean and standard deviation of the experimental frequencies found for each mode i as shown in Table 5.This confidence interval has been calculated in order to compare experimental and numerical results.The purpose of updating the mechanical properties of the numerical model was to ensure that all mode frequencies were included in the confidence interval.Thus, the updated numerical model was significantly similar to experimental results.The experimental frequencies, confidence intervals, initial and updated model frequencies for Modes 3 and 8 are displayed in Table 6.The new mechanical properties are summarized in Table 3.

Discussions
The present study identified eight modes in the range of 0-150 Hz.Six modes are identified in 80% of the tests.The fifth mode seems difficult to exploit given its low repeatability in the identification.This mode is most certainly due to a measurement artifact.Two of these eight modes (54.55 Hz and 126.86 Hz) are in agreement with the numerical modes, and the identification is supported by an analysis of modal deformations.The parameters used in the updated model clearly match both physical and anatomical realities.This proves that an FE model, once updated, is a possible approach to simulate the mechanical behavior of human lower limbs.It should be noted that these two modes are differentiated by the sensor located in the middle of the tibia.The tibia actually contains a vibration node for solicitations around 118 Hz.On the other hand, the thigh does not have a vibration node in the frequency range of interest.
In order to validate the numerical model, this study proposes a methodology that allows describing the modal behavior of the leg with no direct measures.This methodology can come in addition to those used in the literature that focus on modal parameters of a particular bone.One can distinguish two main in vivo methods that are compared in the work of Cornelissen in the late 1980s [35].The Impulse Response (IFR) method involves generating impact with a shock hammer where the leg is relaxed and left hanging.The Bone Resonance Frequency (BRA) method consists of generating a sinusoidal signal over a frequency range with a vibrating shaker.These methods involve the measurement of the forces and their responses on the bone: it is therefore impossible to carry out measurements under real conditions.An OMA offers an interesting alternative, being an application in real conditions where the input forces are not known.Recent developments in OMA make it possible to identify the modal parameters, with possible applications in various sports with stationary impacts or noise.In terms of instrumentation, the proposed setup (featuring four sensors) extracts two modes in the 0-150-Hz range.The sensor in the middle of the tibia is essential to differentiate modes.It might be wise to deploy a fifth sensor in the middle of the thigh to observe a possible modal deformation displaying a vibration node on the femur.However, properly attaching the sensor would be difficult in view of the presence of the rectus femoris.Experimental studies propose a mathematical modeling of the phenomenon and could be integrated directly into the FE model [36].
Concerning the identified frequencies, the frequency of 52.54 Hz is very similar to that found in previous studies [37,38].Both resonance frequencies were determined by using a vibratory excitation at different frequencies between 10 and 90 Hz.The resonance frequency of the ankle was found to lie in the 40-50-Hz range.The modal deformation of the second mode at 126.9 Hz matches the modal deformations of the tibia, which is identified experimentally in the work of Van der Perre [39] and numerically in that of Taylor [40].The latter highlights the torsional modes obtained for higher values (521 Hz).The protocol of the present study unfortunately prevents the measurement of this deformation.
Considering results in a qualitative way, the frequency values lie in wide ranges.The tibia study shows that the bending resonance frequencies are of the order of magnitude 270 Hz and 340 Hz [39].The modal frequencies for the femur are quite controversial.In free-free boundary conditions, Kumar [41] identified modal frequencies at 449, 524 and 1109 Hz, while Gupta [42] determined frequencies of 722 Hz for bending mode and 1582 Hz for double bending.In fixed-fixed conditions, the extracted frequencies are 1211, 1269 and 2818 Hz [41].These studies identified modes on a particular bone with boundary conditions that are completely different from those depicted in the present paper.Our methodology is dependent on the mass of the subject and the whole anatomical environment, which accounts for the lower values and corroborates previous studies on the influence of anatomical parts.Tsuchikane [43] shows from cadaverous specimens that the resonance frequencies are uninfluenced by the skin, increased with the dissection of the muscles and foot and decreased with the dissection of the femur and the fibula.This is confirmed by the work of Cornelissen [35], who adds that the joint has no influence on the modal frequencies of the tibia, but alters modal damping by 10-16%.These results are in agreement with the present work, where the values of frequencies are much lower than the values determined in vitro.In addition, the results of the sensitivity of the parameters show that a variation of 50% in the parameters of soft tissues, including the trunk and hard tissues, leads to a variation of about 10% in the frequencies.
The proposed FE model is a preliminary model that requires many developments, but the values identified experimentally are very similar to the numerical ones, even without updates.Several perspectives are possible.The model does not integrate the fibula, which may actually be necessary: the works of Cornelissen and Tsuchikane show that the fibula makes the lower limb more rigid.Moreover, the work of Tseng [44] proves that the fibula has modal properties close to those of the tibia (173.5, 249.3 and 516.2 Hz): there is no denying that it has a given influence on the dynamic behavior of the limb.Joint cartilage has also not been modeled: its influence, especially on damping, may be investigated in future works.A possible source of error is also the choice of a linear elastic constitutive law for soft tissues.Switching to hybrid models featuring Ogden-type hyperelasticity and exponential anisotropic behaviors [45] may improve the model.Then, it would be interesting to calculate the resonance modes at different knee angles.The work of Munera [8] shows that the resonance frequencies are a function of the knee angle during a squat.The identified frequencies could be located in the excited frequency range 6-8 or 9-20 Hz, matching that of running activities [46].Finally, upgrading to a 3D model would most certainly provide a greater number of modes in the range 0-150 Hz, two-dimensional modeling only allowing extraction of the modes in the mid-lateral plane.
Reading the standards on human exposure to vibrations, the frequency of 52.54 Hz seemed to be the most harmful frequency for lower limbs.Indeed standard ISO 2631-1 [23] defines a frequency weighting filter: the gain at 52 Hz is 0.87, but drops below 0.10 beyond 150 Hz.In the field of public transport and sports, the development of specific materials able to absorb these specific frequencies will decrease the amount of vibrations received by the human lower limbs.For example, a ground made of this material in public transport or in industries could prevent or at least significantly decrease diseases such as arthritis among workers.It could also drastically reduce the tiredness of lower limbs for passengers and make travel far more comfortable.We can note that Taiar [47] determined the vibration responses of anti-fatigue mats to reduce the pain associated with vibratory exposure by considering these standards and defined the optimal geometry and design for risk prevention.

Figure 1 .
Figure 1.(Left) Accelerometers locations and orientations; (right) geometry of the lower limb model in the Modal Analysis software.Points 1, 2, 3, 4 are the location of the accelerometers and the number of node for the numerical model.

Figure 2 .
Figure 2. Numerical model of lower limb (left); mesh detail of hip joint region (right).

Table 2 .
Anthropometric table for a virtual subject 175 cm tall.

Table 3 .
Initial and updated mechanical properties.E: Young modulus, D: Density.

Table 5 .
Results of the modal identification on eleven subjects (Protocol 2).

Table 6 .
Comparison between the confidence interval and the frequencies calculated by the numerical model.